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Abstract 


Reliable  Spaee  Situational  Awareness  (SSA)  is  a  reeognized  requirement  in  the  eurrent 
eongested,  eontested,  and  eompetitive  environment  of  spaee  operations.  A  shortage  of 
available  sensors  and  reliable  data  sources  are  some  current  limiting  factors  for  maintaining 
SSA.  Unfortunately,  cost  constraints  prohibit  drastically  increasing  the  sensor  inventory. 
Alternative  methods  are  sought  to  enhance  current  SSA,  including  utilizing  non-traditional 
data  sources  (external  sensors)  to  perform  basic  SSA  catalog  maintenance  functions. 
Astronomical  data,  for  example,  routinely  collects  serendipitous  satellite  streaks  in  the 
course  of  observing  deep  space;  but  tactics,  techniques,  and  procedures  designed  to 
glean  useful  information  from  those  collects  have  yet  to  be  rigorously  developed.  This 
work  examines  the  feasibility  and  utility  of  performing  ephemeris  positional  updates  for 
a  Resident  Space  Object  (RSO)  catalog  using  metric  data  obtained  from  RSO  streaks 
gathered  by  astronomical  telescopes.  The  focus  of  this  work  is  on  processing  data  from 
three  possible  streak  categories:  streaks  that  only  enter,  only  exit,  or  cross  completely 
through  the  astronomical  image.  Successful  use  of  this  data  will  aid  in  resolving 
uncorrelated  tracks,  space  object  identification,  and  threat  detection.  Incorporation  of 
external  data  sources  will  also  reduce  the  number  of  routine  collects  required  by  existing 
SSA  sensors,  freeing  them  up  for  more  demanding  tasks.  The  results  clearly  demonstrate 
that  accurate  orbital  reconstruction  can  be  performed  using  an  RSO  streak  in  a  distorted 
image,  without  applying  calibration  frames  and  that  partially  bound  streaks  provide  similar 
results  to  traditional  data,  with  a  mean  degradation  of  6.2%  in  right  ascension  and  42.69% 
in  declination.  The  methodology  developed  can  also  be  applied  to  dedicated  SSA  sensors 
to  extract  data  from  serendipitous  streaks  gathered  while  observing  other  RSOs. 
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LEVERAGING  EXTERNAE  SENSOR  DATA  EOR  ENHANCED  SPACE 


SITUATIONAE  AWARENESS 


1.  Introduction 


1.1  Research  Motivation 

Recent  space  activities,  including  the  2007  anti-satellite  (ASAT)  test  by  China  and 
the  2009  collision  between  a  non-operational  Kosmos  satellite  and  an  operational  Iridium 
satellite,  have  increased  global  concern  for  reliable  Space  Situational  Awareness  (SSA) 
[1,  2].  The  current  U.  S.  National  Space  Policy  requires  timely  and  accurate  global  SSA 
[3].  According  to  this  policy  document,  the  SSA  community  should  utilize  information 
from  commercial,  civil,  national  security  sources,  and  foreign  partners  to  detect,  identify, 
and  attribute  actions  in  space  that  violate  long-term  sustainability  of  the  space  environment. 
The  policy  also  calls  for  the  development  and  modernization  of  SSA  capabilities  [3]. 

In  2013,  General  Shelton,  commander  Air  Eorce  Space  Command  (AESPC),  stated 
that  SSA  provides  the  Air  Eorce  the  ability  to  see  and  understand  threats  on  Earth  and  in 
space,  enabling  successful  operations  in  multiple  domains  [4].  The  Joint  Space  Operations 
Center  (JSpOC)  has  stated  that  one  of  their  limiting  factors  for  the  desired  level  of  SSA 
is  a  lack  of  available  sensors  [5].  Although  all  orbit  regimes  are  a  concern,  the  Air 
Eorce  Research  Eaboratory  (AERE)  is  currently  researching  methods  to  characterize  the 
Geosynchronous  (GEO)  regime.  GEO  challenges  include  dim  and  infrequently  observed 
objects,  making  orbital  determination  difficult  and  sustainable  SSA  elusive  [6]. 

Recent  efforts  by  AERE  have  focused  on  developing  an  automated  system  known  as 
Search  and  Determine  Integrated  Environment  (SADIE).  SADIE  is  a  Multi-Hypothesis 
Eilter  (MHE)  designed  to  automate  the  current  time-intensive  process  of  resolving 
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Uncorrelated  Tracks  (UCTs)  and  updating  the  space  catalog  [7].  UCTs  are  defined  as 
observed  tracks  (or  objects  in  space  moving  across  the  Field-of-View  (FOV)  of  a  sensor) 
that  do  not  correspond  to  known  objects.  SADIE  performs  its  mission  by  inputting 
observational  data:  Time,  Elevation,  and  Azimuth  (TEA)  and  range,  from  both  Electro- 
Optical  and  radar  sensors  around  the  world  and  using  MHE  to  resolve  these  tracks  by 
matching  them  to  objects  whose  position  is  currently  unknown.  Unfortunately  SADIE  is 
limited  by  the  data  provided,  requiring  a  minimum  of  four  observations  for  UCT  resolution, 
with  more  observations  providing  a  higher  confidence  in  the  result  [7].  Also  of  note  is  that 
angles-only  solutions,  provided  by  Electro-Optical  (EO)  sensors  and  including  only  right 
ascension,  declination,  and  time,  are  more  accurate  when  geographically  separated  sites 
are  used  to  gather  the  observations  required  [8,  9].  Again,  this  geographic  separation  is  an 
area  with  which  AESPC  struggles,  as  there  are  only  three  optical  sites  used  by  the  Space 
Surveillance  Network  (SSN).  Eocated  on  Diego  Garcia,  at  Socorro,  New  Mexico,  and 
Maui,  Hawaii  these  sites  provide  limited  views  of  the  sky  and  are  tasked  with  all  AESPC 
required  optical  observations  [10]. 

It  is  apparent  that  a  need  for  more  SSA  sensors  exists,  but  these  sensors  don’t  need  to 
be  dedicated  to  the  SSA  mission.  Currently  the  SSN  uses  several  contributing  sensors  to 
provide  additional  data  sources  as  required,  and  the  Defense  Advance  Research  Projects 
Agency  (DARPA)  is  also  exploring  the  idea  of  using  sensors  part  time  as  part  of  the 
Orbital  Outlook  program.  This  outsourcing  demonstrates  that  EO  sensors  designed  for 
other  missions  can  also  be  tasked  to  collect  SSA  data.  This  work  herein  will  explore  the 
idea  of  using  non-dedicated  SSA  sensors  as  an  additional  data  source  for  enhancing  SSA. 

Many  sensor  systems,  such  as  EO  systems  used  by  astronomers,  routinely  record 
satellites  while  performing  their  intended  scientific  purpose.  This  work  explores  the  idea 
of  using  these  serendipitous  collects  to  improve  the  global  SSA  environment.  Unlike 
contributing  sensors,  or  those  used  by  the  Orbital  Outlook  program,  the  sensors  considered 
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in  this  work  referred  to  simply  as  ‘external  sensors,’  will  be  performing  a  mission  other 
than  SSA  and  providing  any  images  including  conveniently  captured  satellites  for  analysis. 
To  achieve  enhanced  SSA  it  is  desirable  for  these  sensors  to  be  able  to  detect  dim  objects 
from  geographically  separated  locations.  These  sensors  should  provide  the  same  TEA  data 
as  current  dedicated  SSA  sensors. 

1.2  Research  Objective 

Several  potential  SSA  resources  exist  outside  the  Air  Force.  Currently  it  is  unknown 
if  data  from  these  resources  can  be  included  in  the  global  SSA  environment.  This  inclusion 
needs  to  be  accomplished  using  a  method  that  ensures  the  data  is  trustworthy,  usable,  and 
will  not  corrupt  the  data  that  has  already  been  validated  by  the  SSN. 

One  potential  resource  is  the  astronomy  community,  as  the  serendipitous  streaks 
collected  during  astronomical  observations  may  be  very  beneficial  to  the  SSA  mission. 
Furthermore,  other  telescope  networks  that  observe  man-made  space  objects  often  referred 
to  as  a  Resident  Space  Object  (RSO),  such  as  the  Exoanalytic^  network.  Orbital  Outlook, 
or  Air  Force  Institute  of  Technology  (AFIT)’s  TeleTrak  system,  can  also  provide  useful 
SSA  data  [11-13].  Any  increase  in  available  resources  will  also  relieve  the  demands  being 
placed  on  the  SSN  providing  not  only  more  data,  but  also  allowing  the  Ground  Based 
Electro-Optical  Deep  Space  Surveillance  (GEODSS)  sensors  to  focus  on  more  demanding 
tasks,  such  as  UCT  resolution.  Space  Object  Identification  (SOI),  or  threat  detection. 

The  objective  of  this  research  is  to  determine  requirements  for  external  sensor  data  to 
be  of  value  for  SSA.  This  work  focuses  on  addressing  the  concerns  of  the  SSA  community 
for  using  sensors  not  controlled  by  the  Department  of  Defense  (DoD)  by  seeking  to 
determine  a  data  calibration  method  that  allows  for  reliable  data  extraction.  Methods  for 
using  the  data  provided  by  the  network  will  be  explored;  the  ability  to  perform  orbital 

'  Exoanalytic  is  a  company  offering  technical  solutions  to  government  and  private  industry  in  the  areas  of 
space,  missile  defense,  and  research  &  development;  Mission  Viejo  CA 
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determination  seems  obvious,  but  the  data  requirements  for  this  application  are  strict.  Every 
observation  collected  will  contain  both  along-track  and  cross-track  data.  Along-track  error 
is  the  largest  error  source  in  the  current  Two  Line  Element  set  (TEE);  documented  in 
Section  2.3.  This  research  will  seek  to  experimentally  demonstrate  the  extraction  of  right 
ascension,  declination,  and  time  data  from  distorted  images  and  provide  it  for  enhanced 
SSA. 

1.2.1  Research  Questions 

This  research  will  be  broken  into  two  phases,  each  phase  addressing  necessary 
questions  to  determine  the  efficiency  of  external  sensors  for  SSA.  The  first  phase,  presented 
in  Chapter  3,  wifi  seek  to  find  the  best  data  caiibration  technique  for  extracting  angie  data 
from  images.  This  objective  wifi  be  accompiished  by  comparing  severai  existing  techniques 
to  determine  the  sensitivity  of  each  process  to  the  distortions  found  in  EO  data.  The  focus 
of  Phase  1  work  wifi  be  to  answer  the  foiiowing  research  questions: 

1 .  What  is  the  ievei  of  accuracy  achieved  by  each  caiibration  method? 

2.  Which  method  best  accounts  for  image  distortion? 

The  work  of  Phase  11,  presented  in  Chapter  4,  examines  the  feasibility  and  military 
utility  of  performing  in-track  and  cross-track  updates  for  an  RSO  catalog  using  metric 
data  obtained  from  streaks  gathered  by  astronomical  telescopes.  Traditional  orbital 
determination  uses  bounded  streaks,  along  with  image  time  and  exposure  duration,  to 
determine  an  accurate  observation  time  and  develop  position  vectors.  Identification  of 
the  center  of  a  streak  is  used  to  provide  precise  timing  of  the  RSO  [8,  14].  The  focus 
of  this  work  is  on  gathering  and  using  three  possible  data  categories:  streaks  that  only 
enter,  only  exit,  or  cross  completely  through  the  astronomical  image.  Streaks  that  only 
enter  or  only  exit  an  image  have  incomplete  timing  data  and  using  them  relies  on  proper 
end  point  detection.  While  streaks  that  pass  completely  through  the  image  provide  no  time 
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data,  requiring  an  accurate  time  estimate  before  they  can  be  used.  The  focus  of  Phase  II 
work  will  be  to  answer  the  following  research  questions: 

1.  Can  observations  of  tracks  only  partially  in  the  image  be  used  to  improve  SSA? 

2.  Can  observations  of  tracks  passing  fully  through  the  image  be  used  to  improve  SSA? 

1.2.2  Assumptions  and  Limitations 

This  research  assumes  that  the  stars  in  Naval  Observatory  Merged  Astrometric  Dataset 
(NOMAD)  catalog  are  a  representative  subset  of  the  U.  S.  Naval  Observatory  Charge- 
Coupled  Device  Astrograph  Catalog  (UCAC)  used  by  DoD.  Both  catalogs  are  developed 
and  maintained  by  the  U.  S.  Naval  Observatory  (USNO),  but  as  the  NOMAD  catalog  is 
much  easier  to  obtain  than  the  UCAC,  NOMAD  will  be  used  as  the  star  catalog  for  this 
research. 

The  work  will  be  limited  to  observing  GEO  satellites  with  TEE  available  in  the  public 
catalog.  Data  will  further  be  limited  to  what  can  be  collected  by  the  TeleTrak  system  at 
AEIT.  Effort  will  be  taken  to  make  these  images  as  representative  of  larger  systems  as 
possible  but  the  EOV  observed  will  be  2-4x  smaller  than  what  would  be  expected  from 
a  larger  survey  telescope  used  by  astronomers.  Einally,  this  research  will  focus  on  using 
images  provided  by  telescopes  in  sidereal  tracking  mode. 

1.3  Dissertation  Overview 

This  document  is  divided  into  five  chapters  and  four  appendices.  Chapter  2  contains  a 
literature  review  of  topics  relevant  to  the  research  that  follows,  divided  into  eight  sections. 
Section  2.1  defines  basic  terminology  used  in  subjects  relevant  to  the  research.  Section 
2.2  covers  current  SSA  methods.  Section  2.3  discusses  various  challenges  faced  by  the 
community.  Section  2.4  covers  alternative  methods  for  determining  position  vectors  used 
in  astrometry.  Sections  2.5  and  2.6  provide  an  overview  of  influential  and  parallel  research 
that  both  benefit  and  enable  aspects  of  this  work,  further  demonstrating  that  the  focus  of  this 
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research  should  be  on  the  validation  of  the  data  produced  by  the  sensors.  Section  2.7  covers 
the  basics  of  orbital  determination.  Section  2.8  discusses  the  usefulness  of  astronomy  data 
and  missions. 

Chapter  3  provides  the  method  and  analysis  used  to  determine  the  best  astrometry 
approach  for  processing  Electro-Optical  data  in  the  experiments  that  follow.  Chapter  4 
covers  the  experiments  used  to  determine  the  usefulness  of  partial  and  full  streaks  in 
imagery  and  provides  the  analysis  of  the  experimental  results.  Chapter  5  reports  the 
findings  of  this  research  and  draws  conclusions  based  on  the  analysis  and  suggested  future 
work. 

The  document  concludes  with  four  appendices.  Appendix  A  sets  up  the  image 
processing  and  calibration  process  and  its  verification.  Appendix  B  provides  the  validation 
and  verification  of  the  astrometry  models  tested  in  Chapter  3.  Appendix  C  discusses  the 
use  of  Goff’s  algorithms  to  perform  Unscented  Kalman  Filter  (UKF)  orbital  determination 
used  to  determine  the  usefulness  of  full  RSO  streaks  [15,  16].  Finally,  Appendix  D  covers 
current  options  for  tools  that  can  be  used  to  calibrate  the  images,  extract  orbital  elements, 
and  return  meaningful  results. 
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2.  Background 


This  chapter  is  designed  to  motivate  the  research,  define  terms,  and  establish  the 
context  of  similar  scholarly  work.  The  chapter  begins  by  providing  important  definitions 
within  the  SSA  community.  An  exploration  of  current  SSA  techniques  will  then  be 
provided  followed  by  a  discussion  of  work  conducted  in  the  past  as  well  as  ongoing  work 
related  to  this  effort.  Finally,  the  chapter  will  explore  alternative  methods  for  conducting 
some  of  the  critical  aspects  in  SSA  that  exist  outside  the  community  but  may  be  useful  for 
image  calibration. 

2.1  Space  Situational  Awareness  Common  Terms 

2.1.1  Space  Situational  Awareness  Defined 

According  to  the  Space  Foundation,  SSA  refers  to  the  ability  to  “view,  understand  and 
predict  the  physical  location  of  natural  and  man  made  objects  in  orbit  around  the  Earth” 
[17].  AFSPC  SSN  Site  Information  Handbook  defines  SSA  as  “the  battle-space  awareness 
required  for  planning,  executing,  and  assessing  protection  of  space  assets,  prevention  of 
hostile  actions,  and  negation  of  hostile  resources  in  all  mediums.”  The  handbook  goes  on  to 
state  that  SSA  is  a  result  of  having  sufficient  knowledge  about  current  and  planned  “space- 
related  conditions,  constraints,  capabilities  and  activities”  [10].  The  National  Space  Policy 
states  that  SSA  is  created  using  “space  surveillance,  intelligence  and  other  information”  [3]. 
Recent  activities  on  orbit,  including  the  2007  ASAT  test  by  China  and  the  2009  collision 
between  a  non-operational  Kosmos  satellite  and  an  operational  Iridium  satellite,  have 
increased  global  concern  for  reliable  SSA.  Space  surveillance  includes  various  methods 
of  detection  and  tracking  of  RSOs  using  astrodynamics  techniques  developed  in  the  last  50 
years. 
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2.1.2  Astrodynamics 

The  field  of  astrodynamics  seeks  to  apply  the  laws  of  celestial  mechanics  to  the 
orbits  of  man-made  objects  around  the  Earth  or  other  celestial  bodies.  These  orbits  are 
fundamentally  derived  using  Newtonian  mechanics  but  also  must  account  for  perturbing 
forces  acting  upon  these  bodies.  Perturbing  forces  include  atmospheric  drag,  solar  radiation 
pressure,  and  gravitational  forces  that  that  require  more  detailed  analysis  than  needed  for 
celestial  objects  [6]. 

2.1.3  J2000  Reference  Epoch 

In  both  celestial  and  orbital  mechanics,  an  epoch  is  a  reference  point  in  time  for 
some  time-varying  quantity,  including  orbital  elements  and  celestial  coordinates  as  they 
are  subject  to  time  varying  perturbations  [18].  The  primary  function  of  an  epoch  is  to 
specify  the  parameters  of  motion.  An  epoch  is  always  referenced  to  a  set  time  such  as  the 
Julian  epoch  (J2000)  the  standard  epoch  in  use  today  by  astronomers.  The  Julian  epoch  is 
exactly  1200  terrestrial  time  on  January  1,  2000  in  the  Gregorian  calendar  [19].  This  epoch 
will  be  used  in  Chapters  3  and  4  to  define  the  motion  of  both  the  stars  and  satellites. 

2.1.4  Earth  Centered  Inertial  Coordinate  System 

Also  important  in  defining  the  location  of  RSOs  is  the  reference  system  used.  The 
preferred  reference  system  for  reporting  satellite  locations  to  the  JSpOC  is  the  Mean 
Equator,  Mean  Equinox  (MEME)  ECI  system.  The  ECI  frame  is  defined  using  the  Earth’s 
orbital  plane  and  the  orientation  of  the  Earth’s  axis  of  rotation  [8],  where  the  first  axis  points 
in  the  direction  of  the  vernal  equinox,  the  second  axis  is  90°  east  along  the  ecliptic  plane, 
and  the  third  axis  is  along  the  Earth’s  axis  of  rotation. 

In  this  coordinate  system  the  position  of  both  satellites  and  stars  are  given  in  right 
ascension  (a),  and  declination  (d).  Right  ascension  is  defined  to  be  measured  positive  to 
the  east  in  the  plane  of  the  equator  from  the  direction  of  the  vernal  equinox;  although  it  can 
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be  reported  from  0^’  -  24^,  in  this  work  degrees,  0°  -  360°,  are  used.  The  declination  is 
measured  northward  from  the  equator  (0®  -  90°)  [8]. 

Throughout  this  work  unless  otherwise  noted,  satellite  positions  will  be  given  in  ECI 
right  ascension  and  declination.  Other  coordinate  systems  used  will  be  discussed  as  they 
appear  throughout  the  work  along  with  the  required  transformation  to  the  ECI  coordinate 
system. 

2.1.5  Uncorrelated  Tracks 

The  term  Uncorrelated  Track  (UCT)  is  a  reference  to  tracked  objects  that  do  not 
correlate  with  known  catalog  objects.  They  are  most  commonly  created  through  on-orbit 
breakups,  maneuvering  satellites,  and  small  objects  that  are  only  occasionally  tracked. 
Each  day  numerous  UCTs  occur  and  performing  correlation  of  these  objects  to  create 
the  required  ephemeris  data  for  catalog  maintenance  is  currently  a  manual  process  [20]. 
Recent  efforts  by  Air  Eorce  Research  Eaboratory  (AERE)  are  attempting  to  automate  this 
process  using  a  program  called  Search  and  Determine  Integrated  Environment  (SADIE), 
see  Section  2.6.2.  Eor  more  on  UCTs  as  described  in  Section  2.3.4. 

2.1.6  Flexible  Image  Transport  System 

The  most  common  type  of  file  produced  by  both  observational  astronomy  and  non- 
SSN  EO  sensors  is  a  Elexible  Image  Transport  System  (fits)  file-format.  The  fits  file  is 
composed  of  a  plain  text  header  followed  by  a  binary  image.  The  header  at  a  minimum 
includes  information  coordinates  of  the  image  source  location  in  World  Coordinate  System 
(WCS).  The  header  may  also  contain:  target,  telescope,  astronomer,  date  observed, 
wavelength  used,  and  any  data  reduction  steps  taken.  The  presence  of  this  data  depends 
on  the  equipment  used  to  collect  the  image  and  the  settings  input  by  the  observer.  The 
binary  image  can  be  composed  of  a  single  image  plane,  or  many  planes  of  images.  The 
‘image’  can  also  be  a  spectrum  or  stack  of  spectra,  or  a  table  of  data  [21]. 
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A  fits  image  is  incredibly  data  rich,  consisting  of  16  to  32  bits  per  pixel.  In  other  words 
there  are  at  least  2^^  (65,536),  and  possibly  as  many  as  2^^  (4.3  billion),  discrete  levels  of 
data  for  each  pixel.  Because  of  the  level  of  information  contained  in  each  pixel,  it  is  worth 
noting  that  saving  a  fits  as  another  file  type  results  in  the  loss  of  considerable  data.  For 
example  a  Joint  Photographic  Expert  Group  (jpg)  file  format  only  stores  2^  (256)  discrete 
levels  of  data.  Therefore  saving  a  fits  file  as  a  jpg  results  in  the  compression  of  data  from 
65,536  levels  to  merely  256  levels  [22].  This  work  will  focus  on  using  fits  images  provided 
by  astronomical  sensors  to  prevent  the  data  loss  that  occurs  when  changing  file  types. 

2.2  Current  Space  Situational  Awareness  Sensors 

Currently  AFSPC  uses  the  SSN  as  a  primary  means  of  achieving  SSA.  The  SSN 
is  comprised  of  phased  array  radars,  mechanical  radars,  optical  systems  and  a  space-based 
optical  system.  In  the  past  the  SSN  also  had  a  multistatic  radar,  but  that  was  shut  down  on  1 
October  2013  [23].  These  systems  are  covered  in  detail  in  the  AFSPC  SSN  Site  Information 
Handbook  and  will  only  be  discussed  in  this  document  as  required  for  understanding  in  the 
context  of  the  research[10].  Other  countries  also  have  systems  that  contribute  to  global 
SSA.  Weeden  et  al.  provides  an  excellent  survey  of  these  systems  that  can  be  found  in 
[24]. 

2.2.1  Optical  Tracking  Techniques 

Optical  tracking  can  be  broken  into  two  methods:  Sidereal  track  and  Rate  track.  The 
Sidereal  track  method  employed  by  telescopes  is  accomplished  when  the  telescopes  FOV 
moves  across  the  sky  at  the  same  rate  as  the  stars  appear  to  move.  This  tracking  method 
keeps  the  star  field  static  in  the  image  taken;  however  any  satellites  in  the  image  will  appear 
as  streaks  of  light  across  the  image.  When  a  GEODSS  telescope  employs  this  method  it 
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can  use  the  fixed  star  image  for  photometric^  calibration,  as  described  in  Section  2. 2. 2.3 

[10]. 

Rate  track  methods  are  performed  by  having  the  telescope  move  to  keep  the  satellite 
stationary  in  the  image.  When  this  method  is  utilized,  the  satellite  appears  as  a  single  white 
dot  and  the  stars  appear  as  streaks  in  the  image  plane.  Both  the  sidereal  and  rate  track 
methods  are  employed  by  GEODSS  to  image  satellites  to  collect  data[10]. 

Traditional  orbital  determination  techniques  use  streaks  bounded  within  the  image, 
along  with  image  time  and  exposure  duration,  to  determine  an  accurate  observation  time 
and  corresponding  position  vectors.  Identification  of  the  midpoint  of  a  streak  is  used  to 
provide  precise  timing  of  the  RSO  as  it  correlates  to  the  middle  of  the  exposure  which  can 
be  determined  with  high  accuracy  [8,  14,  25]. 

GEODSS  are  dedicated  SSN  sensors  tasked  with  providing  timely  astrometric  and 
photometric  data  on  tasked  deep  space^  objects  in  support  of  the  Space  Surveillance 
mission.  Currently  there  are  three  operational  GEODSS  sites  located  in  Socorro,  NM; 
on  the  island  of  Diego  Garcia  in  the  British  Indian  Ocean  Territories;  and  finally  on  the 
island  of  Maui,  HI  [10].  This  distribution  of  telescopes  provides  complete  coverage  of  the 
GEO  belt  [24].  Each  of  these  sites  has  three  telescopes  with  a  40-inch  aperture  and  a  two- 
degree  EOV.  The  optical  system  is  capable  of  detecting  objects  10,000  times  dimmer  than 
the  human  eye  and  is  naturally  limited  by  local  weather  conditions  [10].  These  telescopes 
provide  the  majority  of  all  deep  space  and  GEO  tracks  per  day  for  the  JSpOC  SSA  Cell 
[10]. 

2.2.1. 1  Dual  Angle  Electro-Optical  Tracking 

Both  Gauss  and  Eaplace  developed  methods  for  determining  the  orbit  of  distant 
objects  based  solely  on  visual  observations,  a  serious  challenge  in  mathematical  astronomy 

^The  goal  of  photometric  calibration  is  to  accurately  measure  the  intensity  (per  wavelength)  of  the  light 
captured 

^RSOs  are  considered  to  be  in  deep  space  if  there  orbital  period  is  greater  than  225  minutes  [26] 
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at  the  time.  Today  these  methods  are  widely  used  to  determine  the  orbits  of  satellites.  Dr. 
Pedro  Ramon  Eseobal  provides  detailed  solutions  to  both  methods  in  his  book  Methods  of 
Orbit  Determination  [14]. 

The  basie  eoneept  of  the  Gaussian  method  is  to  take  a  minimum  of  three  line  of 
sight  veetors  and  their  eorresponding  times  and  eonvert  them  to  the  three  position  veetors 
required  to  “fit”  an  orbit.  If  more  than  three  veetors  are  available,  the  solution  ean  be 
improved  using  a  number  of  algorithms.  For  the  eonversion  to  be  aeeurate,  several  site 
eharaeteristies  must  be  known  ineluding  the  observer’s  Loeal  Sidereal  Time  (LST),  latitude 
and  longitude,  and  altitude  as  well  as  the  target  azimuth  and  elevation.  For  ease  of 
ealeulations  all  measurements  are  eonverted  to  the  ECI  eoordinate  system. 

Unfortunately,  no  system  is  perfeet  and  the  data  of  dual  angle  observations  will  not 
be  without  error.  Therefore  it  is  neeessary  to  note  the  level  of  preeision  and  aeeuraey 
provided  by  the  equipment  used  for  these  ealeulations.  Typieally,  as  pointed  out  by  Wiesel, 
systems  are  designed  by  finding  and  eliminating  all  large  error  sourees  until  the  eost  of 
doing  so  exeeeds  the  benefit  of  inereased  aeeuraey  [9].  The  eost  to  aeeuraey  ratio  is  a 
design  parameter  driven  by  requirements;  therefore,  different  teleseopes  will  have  varying 
aeeuraey  that  needs  to  be  aeeounted  for  when  using  an  optieal  system  for  SSA. 

2.2.1.2  Electro-Optical  Error  Sources 

A  Strehl  ratio  is  defined  as  a  measure  of  the  quality  of  optieal  image  formation.  This 
eoneept  was  originally  proposed  by  Karl  Strehl  and  bears  his  name.  The  Strehl  ratio  is 
used  in  situations  where  optieal  resolution  has  been  eompromised  by  any  error  sourees. 
The  ratio  always  has  a  value  between  zero  and  one  with  a  perfeet  optieal  system  reeeiving 
a  one  [27].  Error  sourees  for  optieal  images  inelude  atmospherie  distortions,  instrument 
distortions,  and  sensor  distortions.  Eaeh  of  these  areas  are  diseussed  below. 

The  Earth’s  atmosphere  is  the  largest  eontributor  to  distortion  of  any  system  looking 
into  spaee  from  the  Earth’s  surfaee.  Some  elouds  eompletely  bloek  the  sensor’s  FOV  while 
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others  may  simply  “dim”  the  light  passing  through  them.  Furthermore  water  vapor  and  dust 
also  affect  light  in  time  varying  ways  while  still  allowing  light  to  reach  the  sensor.  Finally 
higher-order  atmospheric  effects  such  as  scintillation,  moment  to  moment  variations  in  the 
atmospheric  refraction,  and  optical  density  changes  in  airflow  past  the  aperture  also  impact 
the  observed  point  spread  function  [27]. 

The  atmosphere  can  be  accounted  for  using  an  atmospheric  model,  such  as 
AFIT’s  Laser  Environmental  Effects  Definition  and  Reference  (EEEDER)  program  [28]. 
Alternatively  if  the  data  is  available,  the  current  atmospheric  effects  can  be  determined 
using  the  stars  available  in  the  image  in  the  same  way  that  Nishimoto  did  in  the  late 
1990’s  [29].  Another  option  is  to  electronically  shift  images  from  one  Charge-Coupled 
Device  (CCD)  cell  to  another  to  counteract  atmospheric  blur,  such  a  system  is  currently 
used  by  Panoramic  Survey  Telescope  and  Rapid  Response  System  (Pan-STARRS)  [30].  To 
illustrate  the  effect  on  possible  astrometry  methods,  the  impact  of  atmospheric  turbulence 
on  an  optical  system  is  modeled  in  Section  3.2.4. 

Sources  of  instrument  distortion  include  jitter,  obscuration  Strehl,  and  higher-order 
aberrations:  defects  in  optics  manufacturing,  and  the  heating/cooling  of  optical  elements 
[27].  Typically  an  optical  system  will  include  measures  to  minimize  some  of  these  effects. 
Other  effects  are  often  accounted  for  by  post-processing  the  data.  The  three  standard  levels 
of  correction  include  raw  uncorrected  data  (Eevel  0),  optical  and  sensor  corrected  (Eevel 
lA/B),  and  ortho-rectified  corrected  where  straight  lines  should  be  straight  and  a  given 
separation  angle  between  two  stars  is  matched  in  the  image  [31].  To  aid  in  the  determination 
of  the  best  astrometry  method  for  this  work  the  impact  of  lens  aberrations  and  jitter  are 
modeled  in  Sections  3.2.2  and  3.2.3,. 

Optical  sensors  used  for  both  SSA  and  astronomy  are  designed  and  maintained  to 
minimize  the  impact  of  the  above  distortion  effects.  Known  as  diffraction- limited  systems 
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these  sensors  have  the  ability  to  produce  images  with  angular  resolution  as  good  as  the 
instrument’s  theoretical  limit  [32]. 

CCD  technology  is  the  preferred  data  generation  method  in  both  astronomy  and  SSA 
due  to  the  CCD’s  sensitivity  and  noise  characteristics.  Unfortunately  these  devices  are  not 
perfectly  uniform;  quantum  efficiency  varies  slightly  from  pixel  to  pixel,  dark  current  is 
generated,  and  only  a  portion  of  each  pixel  is  actually  photo-reactive. 

Photons  arrive  randomly,  during  any  particular  sub  interval,  therefore  the  photons  that 
arrive  during  an  exposure  may  by  greater  or  less  than  the  mean  rate.  This  phenomenon 
is  known  as  Poisson  Noise.  When  dealing  with  faint  objects,  enough  photons  need  to 
be  collected  to  make  a  good  image,  thus  some  finite  integration  time  is  required.  Long 
integration  times  can  obscure  real  variations  in  the  flux  of  photons  from  a  source  impacting 
the  light  curve  ^  data  [33]. 

There  are  also  differences  on  the  order  of  1%  from  one  photo-cell  to  the  next 
between  number  of  incoming  photons  and  the  number  of  electrons  liberated.  This  quantum 
efficiency  effect,  known  as  non-uniformity,  occurs  as  the  photon  signal  is  converted  into 
electrons  and  needs  to  be  accounted  for  [33].  Another  source  of  error  is  the  creation  of  free 
electrons  by  a  CCD  (dark  current),  and  varies  with  temperature  roughly  doubling  every  6° 
Celsius  due  to  thermal  agitation  [33]. 

CCDs  are  adjusted  to  give  an  output  slightly  above  zero  when  no  signal  is  detected. 
Therefore  each  CCD  has  a  zero  point  bias.  Each  pixel  in  the  image  contains  a  measurement 
of  the  system  electrical  offset,  which  is  called  the  signal  bias  [33]. 

Quantization  arises  from  digitizing  the  analog  signal.  This  effect  imposes  a 
“granulation”  on  how  accurately  we  know  the  signal  coming  out  of  the  CCD  camera  [33]. 

The  attachment  of  a  CCD  to  the  optical  system  also  creates  concerns.  Vignetting 
occurs  because  the  optical  system  does  not  illuminate  the  outer  edges  of  the  CCD  as  fully 

^Photometric  measurements  of  an  RSO  are  commonly  referred  to  as  the  ’light  curve’ 


14 


as  it  does  the  center.  Finally  dust  shadowing  caused  by  bits  of  dust  on  a  filter  or  window 
that  cast  shadows  on  the  sensor  affect  the  image  [33]. 

Methods  to  calibrate  sensor  output  exist  for  all  of  these  sensor  effects,  commonly 
referred  to  as  dark  fields,  bias  frames,  and  flat  field  frames.  When  properly  collected,  these 
frames  correct  the  variations  (non-uniformity)  from  the  sensor,  correct  for  the  bias  in  the 
sensor  and  correct  for  any  variations  in  illumination  created  by  the  optical  system  [33]. 

2.2.2  Tracking  by  Determining  Most  Probable  Shape 

Another  application  of  EO  to  the  field  of  orbital  mechanics  is  to  fuse  shape  information 
with  orbital  information  to  create  a  more  accurate  future  position  prediction.  Crassidis  has 
done  considerable  work  in  this  field  and  his  techniques  are  discussed  below. 

2.2.2. 1  Why  does  Shape  Matter? 

Knowledge  of  a  space  object  attributes,  such  as  mass  and  shape,  enables  more  accurate 
non-conservative  force  models.  RSO  attributes  can  also  help  uniquely  identifying  UCT 
or  assist  in  discriminating  between  two  space  objects  on  similar  tracks.  Identification  of 
properly  cataloged  space  objects  is  crucial  to  accurately  predict  future  behavior,  preventing 
further  UCTs,  and  is  critical  for  predicting  conjunctions  [34]. 

At  GEO  the  largest  non-conservative  acceleration  force  is  caused  by  Solar  Radiation 
Pressure  (SRP).  This  force  is  not  guaranteed  to  be  directly  away  from  the  sun  as  size,  shape, 
attitude,  and  material  properties  all  affect  how  satellites  are  affected  by  SRP.  Unfortunately 
these  attributes  may  all  be  currently  unknown,  and  may  change  as  they  age  and  re-orient 
requiring  them  to  be  determined  as  a  function  of  time  [34] . 

2.2.2.2  Attribute  Estimation 

According  to  Crassidis  the  objective  of  optical  SOI  is  to  “identify,  attribute,  and 
track  characteristics  of  space  objects  from  photometric  data.”  This  method  allows  for  the 
characterization  of  RSO  to  determine  the  properties  of  the:  shape,  surface,  and  material 
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along  with  dynamic  states.  Dynamic  states  can  be  used  to  classify  RSOs  further  as  either 
uncontrolled  (debris  or  rocket  body)  or  controlled  (active  satellite)  [34] . 

2.2.2. 3  Methods  to  Determine  Shape 

Numerous  estimators  exist  for  estimating  the  shape  of  an  observed  object.  These 
methods  include  radar  cross-sectioning  approaches,  as  discussed  by  Toru  Sato  et  al.  in 
[35],  and  range  Doppler  interferometry,  which  is  discussed  in  detail  in  [36].  In  Low  Earth 
Orbit  (LEO)  some  objects  can  even  be  imaged  by  ground  based  telescopes.  Often  theses 
methods  are  limited  by  the  size  and  distance  of  the  RSO  to  the  observer  [37]. 

Laser  Radar  (LADAR),  which  returns  a  three  dimensional  scan  of  an  object,  has  also 
successfully  been  used  to  return  the  shape  of  an  object.  An  advantage  of  using  LADAR 
is  the  detail  of  the  shape  returned.  Unfortunately  the  range  of  LADAR  limits  its  use  to 
proximity  or  rendezvous  operations  [38]. 

Lichter  and  Dubowsky  cover  an  application  for  LADAR  proximity  operations  for 
shape  matching  that  is  of  interest  here  both  for  its  use  of  kinematic  data  fusion  and  because 
it  determines  orbital  characteristics  as  part  of  the  solution  method.  Lichter  and  Dubowsky 
take  range  image  data  from  a  synchronized  team  of  LADAR  or  stereo  cameras  and  reduce 
it  into  coarse  kinematic  information  at  each  sample  time.  This  data  is  then  passed  to  a 
Kalman  filter,  which  forms  the  core  of  the  estimation  architecture  [39] . 

The  Kalman  filter  extracts  the  full  dynamic  state  and  inertial  parameters  of  the  object 
using  the  surrogate  measurements,  along  with  an  accurate  dynamic  model.  The  Kalman 
filter  provide  estimates  of  both  rotational  and  translational  positions  and  velocities  of 
the  RSO.  Once  accurate  information  about  the  trajectory  of  the  target  is  known,  shape 
estimation  reduces  to  a  classic  stochastic  map-building  problem  [39]. 

As  the  target  motion  has  been  determined  with  respect  to  the  sensors,  pixel-level  data 
can  be  fused  into  a  probabilistic  map  of  the  target’s  shape  using  sensor  noise  models.  This 
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method  can  be  used  for  small  satellites,  debris,  and  large  flexible  space  structures,  but  as 
mentioned  above,  the  LADAR  sensor  has  a  significantly  limited  range  [39] . 

An  alternative  method  to  LADAR  which  uses  natural  illumination  is  to  observe 
and  measure  the  reflected  sunlight  from  an  RSO.  Sunlight  reflected  off  of  a  satellite  is 
fundamentally  different  than  direct  sunlight.  The  fundamental  difference  in  the  wavelength 
of  reflected  sunlight  can  be  attributed  to  the  reflectivity  of  the  surface  materials,  the  object’s 
size  and  orientation,  and  the  medium  the  light  travels  through  along  the  line  of  sight  to  the 
observer.  This  fact  can  be  leveraged  to  determine  one  or  more  of  these  properties  from  the 
light  curve  of  an  observed  object  [38].  In  fact,  astronomers  have  been  using  light  curve 
data  to  study  asteroid  characteristics  for  some  time  [40]. 

It  is  crucial  when  working  with  light  curve  data  that  the  observer  accounts  for 
environmental  factors  that  impact  the  light  curve.  Most  notably  for  RSOs,  the  atmosphere 
must  be  accounted  for.  In  a  study  performed  by  Nishimoto  et  al.  from  1999  to  2001, 
a  method  was  developed  to  account  for  the  atmospheric  impact  on  light  curve  data. 
Background  measurements  of  well  understood  visible  stars  where  taken  each  night  of 
observations;  these  observations  allowed  the  system  to  account  for  atmospheric  and 
instrument  affects  on  the  reflected  light  [29].  Once  the  external  affects  were  accounted 
for  and  removed  the  wavelength  of  the  reflected  light  was  known.  The  reflected  light 
can  be  compared  to  the  absorptivity  of  known  satellite  materials  to  determine  the  material 
surface  that  reflected  the  sunlight  as  was  demonstrated  by  Cauquy  et  al.  in  2006.  These 
measurements  can  also  be  used  to  identify  two  closely  spaced  spacecraft  built  out  of 
different  material  (aluminum  and  gold  solar  blankets  for  example)  [41]. 

Linares  et  al.  presented  a  method  in  2013  using  an  Multiple-Model  Adaptive 
Estimation  (MMAE)  approach  to  identify  the  most  probable  shape  of  an  orbiting  object 
by  comparing  it  to  a  catalog  of  candidate  shape  models.  MMAE  uses  a  set  of  parallel 
filters,  each  hypothesizing  its  own  shape,  to  determine  the  probable  shape  of  the  object  in 
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question.  The  filters  used  compared  a  known  shape  to  light  curve  data  collected  by  an  EO 
system,  the  details  of  which  are  given  in  Reference  38. 

The  coupling  between  RSO  attitude  and  non-conservative  accelerations,  SRP  and 
atmospheric  drag,  can  be  exploited  to  aid  in  the  estimation  of  an  object’s  trajectory. 
SRP  acceleration  is  a  function  of  the  same  properties  provided  by  light  curve  methods, 
such  as  the  one  developed  by  Linares  et  al.  Therefore  these  methods  also  allow  for  the 
determination  of  the  object’s  inertial  orientation  and  trajectory.  While  the  parallel  filters 
are  checking  the  shape  they  also  simultaneously  reduce  the  passively  collected  EO  data 
to  infer  the  unknown  state  vectors  composed  of  the  object’s  position,  inertial-to-body 
orientation,  and  respective  temporal  rates  [38].  The  Linares  et  al.  method  is  preferred 
to  those  mentioned  above  because  it  can  be  used  on  dim  objects  in  GEO,  relying  only  on 
light  curve  data. 

As  covered  above  several  optical  techniques  are  currently  in  use  in  the  SSA  field. 
Ground  and  space  based  sensors  are  used  to  track  satellites  both  through  angles  only 
techniques  and  shape  estimation.  These  techniques  also  face  challenges  ranging  from  the 
accuracy  of  orbital  positions  to  infrequently  observed  dim  objects. 

2.3  Current  Space  Situational  Awareness  Challenges 

The  National  Research  Council  (NRC)  recently  found  that  current  software  suites  used 
to  perform  astrodynamics  calculations  cannot  meet  current  requirements  to  perform  critical 
tasks  such  as  conjunction  analysis,  and  furthermore,  they  struggle  with  even  fundamental 
tasks  such  as  statistical  orbit  determination.  The  current  system  also  makes  it  difficult  to 
incorporate  changes,  including  those  that  would  enhance  the  system.  The  system  in  use  by 
JSpOC  is  also  not  flexible  enough  to  meet  all  user  needs  [6].  This  section  summarizes  the 
findings  of  the  NRC  in  their  report  titled  ‘Continuing  Kepler’s  Quest-Assessing  Air  Eorce 
Space  Command’s  Astrodynamics  Standards,’  the  full  report  can  be  found  in  Reference  6. 
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As  future  sensors  are  able  to  resolve  smaller  objects  the  problem  of  correlating  each 
new  object  will  increase  leading  to  a  growing  number  of  UCTs.  Many  commercial 
companies  have  demonstrated  innovative  algorithms  that  address  specific  issues  faced  by 
the  SSA  community.  Unfortunately  none  of  the  alternatives  addresses  the  full  requirement 
set  [6]. 

2.3.1  Current  Uncertainty  in  Low  Earth  Orbit  Covariance 

According  to  Osweiler,  using  just  the  publicly  available  TLE  generated  by  the  current 
SSN,  LEO  satellites  have  their  largest  prediction  error  in  the  in-track  velocity  component 
of  the  satellite  position  vector.  This  error  ranges  from  130  km^  to  4600  kirp-.  Osweiler 
attributes  this  error,  both  magnitude  and  direction,  to  atmospheric  drag  [42]. 

2.3.2  Current  Uncertainty  in  Geosynchronous  Covariance 

Accuracy  of  TLE  data  for  Geostationary  orbits  and  specifically  Intelsat  satellites  have 
been  studied  in  great  detail.  Chan  and  Navarro  compared  the  publicly  accessible  TLE  with 
INTELSAT  proprietary  orbital  elements  determined  using  their  own  dual  station  ranging 
system.  The  purpose  of  their  comparison  was  to  find  out  how  reliable  the  TLE  was  for 
conjunction  detection  analysis.  After  comparing  data  for  a  full  year  Chan  found  that  the 
orbital  differences  showed  an  average  of  30-60  km  difference  in  position,  when  the  most 
recent  TLE  was  propagated  to  the  observational  epoch  [43].  Eurther  studies  have  been 
performed  on  the  accuracies  of  GEO  TLEs  in  general.  Eriih  and  Schildknecht  performed  a 
four  year  study  and  after  evaluating  13  GEO  objects  found  an  average  of  25  km  along-track 
and  10  km  cross-track  difference  after  propagating  the  most  current  TLE  to  the  current 
epoch.  They  also  found  that  the  error  grows  the  further  from  the  TLE  epoch,  reaching 
hundreds  of  kilometers  of  track  error  in  just  a  couple  of  days  [44]. 

Osweiler’s  study  shows  that  similar  to  LEO  satellites,  GEO  satellites’  largest  observed 
prediction  error  from  the  TLE  is  also  in-track  or  along  the  velocity  component,  compared  to 
LEO  satellites,  Osweiler  found  that  GEO  satellites  have  a  much  smaller  covariance  matrix 
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overall.  The  magnitude  of  error  varies  by  orbit  as  well.  Satellites  in  stable  orbits  with 
negligible  external  forces  have  an  average  covariance  uncertainty  of  approximately  1  kni^. 
While  satellites  affected  by  gravity  variations  or  SRP  have  a  covariance  range  from  less 
than  1  krn^  to  greater  than  25  krn^  [42].  Recent  work  performed  by  Sabol  et  al.  for  AFRL 
found  similar  results  in  covariance  uncertainty  [7,  45]. 

2.3.3  High  Area-to-Mass  Ratio  Objects 

An  important  category  of  GEO  objects  is  the  High  Area-to-Mass  Ratio  (HAMR) 
object.  HAMR  objects  tend  to  have  area-to-mass  ratios  between  1  and  40  tr?  j kg  [46]. 
Schildknecht’s  study  on  HAMR  objects  revealed  that  the  inclination  and  eccentricities  are 
both  highly  effected  by  SRP,  and  to  a  lesser  extent  solar  and  lunar  gravitational  forces,  with 
eccentricities  ranging  from  0.05  to  0.6  [46].  Schildknecht’s  study  also  showed  that  HAMR 
objects  also  tend  to  have  a  mean  motion  of  about  one  revolution  per  day  about  their  center 
of  gravity  [46] . 

To  acquire  more  information  on  sizes,  shapes,  and  possibly  materials  of  these  HAMR 
objects,  light  curve  data  was  gathered.  The  light  curves  revealed  a  high  variability  of 
information  indicating  both  periodic  and  random  variations  of  several  orders  of  magnitude 
over  a  few  minutes  of  observations  to  constant  brightness  over  durations  greater  than 
10  minutes.  The  data  also  demonstrated  that  the  behavior  of  any  single  object  could 
completely  change  from  one  observation  to  the  next.  These  results  indicate  randomly 
tumbling  objects  with  complicated  shapes  [46]. 

Finally  most  HAMR  objects  seem  to  have  a  relatively  constant  area-to-mass  ratio  over 
time;  however,  some  variations  in  this  ratio  were  recorded  and  could  change  by  up  to  a 
factor  of  four  over  a  few  days  [46]. 

Schildknecht  concluded  that  additional  observing  techniques  would  be  beneficial 
to  determine  the  material  and  the  nature  of  HAMR  objects,  suggesting  photometry, 
spectroscopy,  and  radar  as  potential  methods.  He  also  indicated  that  continuous  tracking  of 
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a  subset  of  high-area-to-mass  ratio  objects  would  be  mandatory  to  ensure  the  space  object 
catalog  was  current  [46] . 

2.3.4  Resolving  Uncorrelated  Tracks 

Coffey  et  al.  make  it  clear  that  correlating  tracks  of  known  space  objects  that  are 
not  associated  with  an  entry  in  the  space  catalog  is  a  complex  and  difficult  problem 
[47].  Currently  the  DoD  catalogs  more  than  22,000  orbiting  objects,  including  active  and 
inactive  satellites,  rocket  bodies,  and  debris  [48].  These  numbers  are  expected  to  grow 
exponentially  when  the  SSN  is  upgraded  in  the  future  allowing  the  detection  and  tracking 
of  even  smaller  objects.  These  smaller  objects  are  already  on  orbit,  but  undetected  by  the 
current  sensors.  When  future  sensors  locate  them  they  will  start  as  UCTs  that  need  to  be 
resolved  compounding  the  current  problem  [7].  A  recent  report  from  the  NRC  found  that: 

“Advanced  data  association  methods  such  as  multiple  hypothesis  tracking  will 
produce  a  significant  improvement  in  the  automation  of  the  space  surveillance 
system,  including  the  sensor  processing.  The  algorithms  should  be  adapted  to 
high-performance  computing  and  advanced  computer  architectures  and  should 
make  use  of  kinematic,  feature,  and  nontraditional  data”  [6]. 

Many  innovative  methods  attempting  to  resolve  UCTs  depend  on  sensor  data 
association  techniques.  Section  2.6.3  discusses  UCT  resolution  techniques  and  Section 

2.3.5  discusses  data  association  methods. 

2.3.5  Data  Association  and  Information  Fusion 

Data  association  and  fusion  of  information  are  fundamental  problems  in  tracking 
multiple  space  objects.  Data  association  becomes  a  problem  when  a  collection  of  reports, 
data  measurements  and  sensor  tracks,  need  to  analyzed  to  determine  which  reports 
originated  from  the  same  object  and  which  belong  to  uncorrelated  tracks  or  false  reports. 
Once  this  correlation  has  been  accomplished,  information  about  a  single  object  from  one 
or  more  sources  needs  to  be  combined  to  improve  our  understanding  of  the  object  [6]. 
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Current  algorithms  are  no  longer  sufficient  to  handle  the  number  of  UCTs,  and  future 
demands  will  far  exceed  current  demand.  Furthermore,  future  demands  will  require  new 
data  association  and  nonlinear  estimation  algorithms,  advanced  computer  architectures, 
orbit  propagation,  and  orbit  determination. 

Currently  AFSPC  uses  a  non-statistical  fixed  gated  data  association  method  known 
as  the  Report  Observation  Association  (ROTAS).  The  ROTAS  system  works  well  with 
widely  spaced  objects  against  a  benign  background.  Recent  developments  have  made  space 
more  congested  with  LEO  breakups  and  clusters  of  objects  at  GEO  that  require  new  data 
association  methods  because  the  objects  can  be  closely  spaced  and  the  association  is  often 
ambiguous  [6]. 

According  to  the  NRC,  a  good  space  object  association  process  should  make  use  of 
other  pieces  of  data  as  well:  Radar  cross  section,  intensity,  and  other  features  could  be  used 
to  identify  an  object  which  will  help  correlate  it  to  an  orbit.  Data  not  traditionally  used  for 
tracking  space  objects  (multi-band  photometry,  radiometry,  etc.),  should  also  be  considered 
[6]. 

A  UCT  poses  an  even  greater  challenge.  The  problem  of  associating  a  UCT  from 
one  sensor  to  another  sensor,  or  even  the  same  sensor  after  several  orbits,  is  quite  difficult. 
This  challenge  defines  the  data  association  problem  today,  considered  a  central  problem  in 
multiple-target  tracking.  The  problem  is  further  complicated  by  miss-association  or  cross¬ 
tagging  of  objects  [6,  49-51]. 

The  data  association  problem  for  space  surveillance  stems  from  the  need  to  partition 
sensor  tracks  into  valid  objects,  and  false  reports  or  UCTs.  Algorithms  for  solving  data 
association  problems  are  generally  divided  into  two  classes:  single-frame  and  multiple- 
frame  methods  [6]. 
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23.5.1  Data  Association:  Single-Frame  Methods 
Single  frame  methods  inelude  Nearest  Neighbor  (NN)  and  Global  Nearest  Neighbor 
(GNN).  Both  NN  and  GNN  work  well  in  light  clutter  and  for  reasonably  spaced  objects^. 
The  method  currently  used  by  AFSPC  is  GNN.  GNN  is  a  two  dimensional  assignment 
problem  which  attempts  to  relate  all  objects  in  a  sensor  fame  to  known  objects  expected  to 
be  present.  Unfortunately  in  the  current  crowded  space  environment  single-frame  methods 
fail  to  appropriately  identify  and  tag  RSOs  leading  to  an  increase  in  UCTs  [6]. 

2.3.5. 2  Data  Association:  Multi-frame  Methods 
Multiple-Hypothesis  Testing  (MHT)  mitigates  miss-association  or  cross-tagging  by 
holding  difficult  association  decisions  in  abeyance  until  additional  information  is  available, 
effectively  improving  current  decisions  by  changing  past  decisions  [6]. 

Although  MHT  methods  have  been  around  for  the  past  30  years,  they  have  only 
started  to  be  used  in  observation  systems  in  the  last  decade.  The  development  of 
similar  advanced  association  methods  for  SSA  will  potentially  address  several  challenges. 
MHT  methods  can  be  adapted  to  advanced  computer  architectures,  with  a  focus  on 
control  of  the  computational  costs  for  these  architectures.  All  available  data  including 
kinematic,  feature,  and  nontraditional  data  (see  Section  2.3.5)  is  required  to  disambiguate 
the  association  process  when  needed.  Association  methods  themselves  should  be  adapted 
to  the  complexity  of  the  problem.  For  widely  spaced  objects  simple  methods  such  as  NN 
can  be  used,  whereas  for  LEO  breakups  and  GEO  clusters,  a  full  MHT  algorithm  may  be 
required  [6] . 

Blackman’s  survey  of  MHT  provides  an  excellent  overview  of  the  basics  of  MHT. 
As  visual  representation  of  the  concept  of  MHT  tends  to  be  the  best  way  to  explain  it, 
Eigure  2.1  is  used  to  show  the  predicted  locations  of  two  targets  and  three  observations 
made  when  the  targets  were  predicted  to  be  present.  With  MHT  any  two  pairs  could  be  the 

^Reasonably  spaced  objects  at  GEO  are  typically  separated  by  at  least  1  km 
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correct  targets,  for  example  0\  and  02  represent  a  single  hypothesis,  called  H\.  Once  the 
observations  are  made  there  are  10  possible  hypotheses.  These  hypotheses  range  from  the 
most  likely  scenario,  P2  belonging  to  02  and  0\  and  03  being  attributed  to  P\,  to  the  least 
likely,  all  three  observations  are  false  alarms  [52]. 


01, 02,  03  s  Observation  positions 
PI,  P2  5  Predicted  target  positions 


Figure  2.1:  Typical  data  association  conflict:  prediction  variance  vs.  observed  objects 


After  Singer,  Sea,  and  Housewright  introduced  the  concept,  Reid  developed  his 
original  algorithm  which  is  implemented  as  a  sequential  k-best,  two  dimensional 
assignment  algorithm.  Reid  defined  a  systematic  approach  for  evaluating  multiple  data 
association  hypotheses  in  a  clutter  or  false  alarm  background  [49,  50,  52-55]. 

Blackman  demonstrates  how  Ried’s  algorithm  works  using  Figure  2.1.  Start  by 
defining  tracks  T\  and  T2  as  the  predicted  positions  P\  and  P2.  Then,  define  all 
newly  formed  tracks,  such  as  T3{P\,OV)  from  the  association  of  FI  with  01,  for  other 
observation  and  track  combinations  to  ensure  all  hypotheses  are  accounted  for.  Finally, 
Reid  also  defines  NT\,  NT 2  and  NT 3  as  the  new  tracks  initiated  from  01,  02,  and  03. 
Now  we  can  define  all  10  hypotheses: 
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Hi  :  n,  T2,  NTl,  NT2,  NTS 

H2  :  T3,  TA,  NTS 

:  TS,  T5,  NT2 


(2.1) 


Compatible  tracks  are  defined  as  those  with  no  common  observations  and  MHT 
hypotheses  are  made  up  of  sets  of  compatible  tracks.  Reid’s  algorithm  approach  also  carries 
over  hypotheses  from  previous  scans.  When  a  new  scan  is  performed,  the  data  are  combined 
and  consolidated  into  only  the  hypotheses  that  meet  the  compatibility  constraint  [52]. 

According  to  Blackman  a  track-oriented  MHT  approach  recomputes  all  hypotheses 
using  newly  updated  tracks  after  every  scan.  Track-oriented  MHT  then  discards  the 
hypotheses  formed  on  the  prior  scan  after  performing  the  update  to  create  new  hypotheses. 
This  method  prevents  the  loss  of  data,  as  the  tracks  that  are  maintained  contain  all  the 
relevant  statistical  data,  but  does  not  maintain  prior  data  beyond  the  current  scan.  Blackman 
prefers  this  approach  when  there  are  many  track  observations  (clustered  environment) 
because  unlike  other  methods,  it  prevents  hypotheses  from  growing  exponentially. 
Alternative  methods  surveyed  by  Blackman  can  be  found  in  Reference  52. 

2.4  Astrometry  Techniques 

The  science  of  astrometry  aims  to  accurately  define  positions  of  objects  in  reference 
to  the  celestial  sphere  [33].  To  accomplish  these  measurements  in  CCD  image  plate 
coordinates,  pixel  locations  corrected  for  distortion  must  be  determined.  If  sky  coordinates 
are  available  for  known  stars  in  the  image  they  can  be  used  to  determine  the  plate 
coordinates.  These  plate  constants  can  then  be  used  to  estimate  the  location  of  objects 
with  unknown  sky  coordinates  [33].  Alternatively  if  camera  parameters  are  known  then  a 
camera  model  may  be  used  to  determine  the  real-world  position. 
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Astrometry  is  practiced  by  astronomers  to  define  the  location  of  newly  discovered 
objects  so  that  they  can  be  accurately  located  by  others.  CCD  images  can  be  utilized  to 
obtain  precise  sky  coordinates  of  celestial  objects.  Other  benefits  of  this  method  include 
the  ability  to  accurately  determine  the  focal  length  of  the  telescope  used  for  taking  the 
image,  as  well  as  the  images  precise  scale  and  orientation  [33]. 

The  use  of  astrometry  to  find  moving  objects  is  also  well  understood.  Techniques  have 
been  developed  for  finding  asteroids  and  comets  and  then  computing  their  orbit  around  the 
sun  [33].  For  man  made  satellites,  the  same  dim  object  detection  techniques  can  be  used  to 
locate  Earth  orbiting  satellites. 

Chapter  3  focuses  on  the  comparison  of  three  methods  of  astrometry  for  angle 
data  extraction:  affine  transform  astrometry,  the  pinhole  camera  model,  and  the  intrinsic 
parameters  camera  model,  also  known  as  the  intrinsic  parameters  calibration  model.  Each 
method  has  its  own  requirements  outside  of  the  image,  and  will  be  discussed  in  detail 
in  Sections  2.4.1  through  2.4.3.  Both  the  pinhole  camera  model  and  the  affine  transform 
method  have  been  successfully  used  in  astrometry  in  the  past  to  track  RSOs  [13,  33].  The 
validation  and  verification  of  each  model  is  covered  in  Appendix  B. 

2.4.1  Affine  Transform  Astrometry 

The  first  method  considered  for  extracting  satellite  position  information  from  an  image 
will  be  the  application  of  an  affine  transform.  Requiring  only  the  proper  identification  of 
stars  in  the  image,  and  number  of  pixels  in  the  frame,  affine  transform  astrometry  requires 
the  smallest  amount  of  data  of  any  of  the  considered  methods  [33]. 

To  use  the  affine  transform  method,  start  by  measuring  {x' ,  y')  locations  for  each 
reference  star  and  collect  the  corresponding  look  angles;  right  ascension  and  declination. 
Then  compute  plate  constants,  transformation  parameters  between  image  pixel  coordinates 
and  sky  coordinates  (a  and  d),  for  the  image.  Once  plate  measurements  are  known,  measure 
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the  {x' ,  y')  position  of  the  target  objects.  Finally,  compute  the  coordinates  of  the  target 
object.  Each  of  these  steps  is  outlined  below. 

2.4.1. 1  Standard  Coordinates 

The  standard  coordinates  used  in  astrometry  assume  that  an  image  of  the  sky  contains 
a  small  enough  FOV  that  the  spherical  surface  can  be  treated  as  a  plane.  Thus,  the  standard 
coordinate  plane  is  defined  as  a  plane  tangent  to  the  celestial  sphere  at  a  point  (ao,  do)  on 
the  sky.  According  to  Berry,  “the  X-axis  is  aligned  with  right  ascension  (a),  the  Y-axis  is 
aligned  with  declination  (d),  and  the  origin  lies  at  the  point  of  tangency.”  Therefore  it  is 
straightforward  to  compute  the  location  (X,  Y)  of  any  star  with  a  known  right  ascension 
and  declination  as  demonstrated  in  Equations  2.2  and  2.3  [33]: 


cos(d)  sin  (a  -  ao) 

cos(do)  cos(d)  cos  (a  -  ao)  +  sin(do)  sin(d) 


(2.2) 


_  sin(do)  cos(d)  cos  (a  -  ap)  -  cos(do)  sin(d) 
cos(do)  cos(d)  cos  (a  -  ao)  +  sin(do)  sin(d) 

These  standard  coordinates  are  useful  as  they  mimic  the  formation  of  images  on  the 
CCD  imager  [33].  Thus  the  left  hand  side  of  Equations  2.2  and  2.3  are  the  same  as  the  left 
hand  side  of  Equations  2.4  and  2.5  allowing  for  the  determination  of  the  plate  constants. 

2.4.1.2  Finding  Plate  Constants 

In  an  ideal  image  the  X  and  Y  axes  used  in  measuring  the  image  would  be  perfectly 
aligned  with  its  right  ascension  and  declination  axes.  The  exact  center  of  the  image  {ap,  do) 
would  be  given.  Under  ordinary  circumstances,  however,  (ao,  do)  will  not  be  perfectly  at 
the  center  and  the  X  and  Y  axes  of  the  detector  will  be  rotated  through  some  angle,  as  seen 
in  Figure  2.2.  These  offsets  will  be  automatically  compensated  for  when  the  mathematical 
transform  from  standard  to  plate  coordinates  occurs  [33]. 
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Figure  2.2:  Plate  eoordinates  vs.  standard  eoordinates 


Offsets  between  the  axes  can  be  caused  by  displacement,  rotation,  tilt,  and  inaccurate 
pixel  dimensions  at  the  CCD’s  operating  temperature.  Instead  of  trying  to  account  for  every 
possible  optic  and  sensor  combination  in  use.  Berry  states  that  it  is  possible  to  write  the 
relationship  between  two  coordinates  systems  as  a  general  affine  transformation  as  shown 
in  Equations  2.4  and  2.5  [33]: 


X  =  ax  +by'  +  c  (2.4) 

Y  =  dx  +  ey  +  f  (2.5) 

With  X,  Y,  x',  y'  as  previously  defined  and  terms,  a-f  the  ‘plate  constants’.  The  plate 
constants  can  be  determined  empirically  from  the  CCD  image  itself  by  measuring  the  x' ,  y' 
location  of  a  minimum  of  three  reference  stars,  and  solving  the  resulting  linear  equations 
for  the  six  unknowns  using  six  independent  equations  as  shown  in  Equation  2.6  [33]: 
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Xi  x[  1  0  0  0  a 

X2  x'2  I  0  0  0  b 

X3  1  0  0  0  c 

=  '  '  (2.6) 

Yi  0  0  0  x[  y[  1  d 

Y2  0  0  0  1  e 

73]  [0  0  0  y'  lj[/ 

If  more  than  three  reference  stars  are  available  in  the  image  they  can  be  used  to 
improve  the  initial  calculation  of  the  plate  constants  by  performing  a  least-squares  solution 
to  produce  the  most  probable  solution  to  an  overdetermined  set  of  equations  [33]. 

2. 4.1. 3  Extracting  Position  Coordinates  from  Streaks 
Once  the  plate  constants  are  known,  the  position  of  the  object  in  the  image  can  be 
determined  from  its  {x' ,  y')  plate  coordinates  utilizing  Equations  2.4  and  2.5  [33].  Finally 
to  convert  the  standard  coordinates  of  the  object  back  to  right  ascension  and  declination  use 
Equations  2.7  and  2.8^  [33]: 


5  =  arcsin 


sin(do)  +  Y  cos(do) 

Vl  +  +  72 


(2.7) 


a  =  ao  +  arctan ( — , .  .  J  (2.8) 

\cos(do)  -t  7  sm(do)/ 

2.4.2  Frontal  Pinhole  Camera  Model 

An  alternative  astrometry  method  is  to  apply  a  pinhole  camera  model  to  the  provided 
image  [13].  This  model  assumes  that  all  light  rays  are  forced  to  go  through  the  optical 
center  which  minimizes  optical  distortions.  A  well-focused  system  designed  such  that 
distortions  in  the  image  are  negligible  effectively  becomes  a  pinhole  camera  [56]. 

For  the  purpose  of  the  comparison  performed  here  the  image  plane  will  be  placed  at 
a  location  equal  to  the  distance  of  the  focal  length  in  front  of  the  origin  (optical  center), 

®The  two  quadrant  inverse  tangent  is  used  in  Equation  (2.8)  and  the  result  is  added  to  ao  to  provide  a. 
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along  the  Z  axis,  as  shown  in  Figure  2.3.  This  plaeement  defines  the  eonventional  frontal 
pinhole  eamera  model  [56]. 


Figure  2.3:  Frontal  pinhole  imaging  model 


When  the  pinhole  eamera  model  is  used,  the  ideal  perspeetive  projeetion  ean  relate 
the  image  of  an  objeet  to  its  physieal  eoordinates  as  shown  in  Figure  2.3  [56].  The  pinhole 
eamera  model,  provided  in  Equation  2.9,  ean  be  used  to  eonvert  the  pixel  loeation  {x' ,  y') 
for  the  observed  RSO  to  an  \XY  ZY  position  veetor. 


x' 

r  0  0, 

A 

= 

0  r  Oy 

Y 

1 

0  0  1 

Z 

(2.9) 


Where  A  in  Equation  2.9  is  a  normalizing  seale  faetor.  As  Figure  2.3  shows  the  sensor’s 
[Z  Y  ZY  plane’s  origin  is  along  the  optieal  axis,  marked  by  the  x  and  y  axis  origin,  plaeing 
it  in  the  middle  of  the  image.  The  \X  Y  ZY  sensor  frame  X  and  Y  eomponents  are  just 
shifted  versions  of  the  pixel  rows  and  eolumns  in  pixel  units,  as  provided  by  Equation  2.10, 
where  N  and  M  are  the  number  of  horizontal  and  vertieal  pixels  of  the  CCD. 


0, 

V+1 

2 

f 

M+1 

L  2  ] 

(2.10) 
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The  downrange  vector  component  (Z)  is  determined  by  the  sensor’s  focal  length,  in 
pixels,  as  derived  from  the  relationship  shown  in  Figure  2.4  [13].  In  Figure  2.4,  6  is  the 
angle  at  which  light  from  the  edge  of  the  FOV  passes  through  the  origin  of  the  coordinate 
system.  The  prime  symbol  is  used  to  indicate  virtual  measurements  inside  an  optical 
instrument  [13,  57].  From  Figure  2.4  you  can  see  that  9'  equals  half  the  angle  that  defines 
the  FOV.  Therefore  the  focal  length  can  be  calculated  in  terms  of  pixels  using  Equation 
2.11^  [13,56]. 


Figure  2.4:  Virtual  focal  length  (/')  as  derived  from  field  of  view  (6) 


N 

= - Tfo^ 

2  arctan 

All  of  the  variables  required  for  Equation  2.9  have  now  been  fully  defined  from  pixel 
coordinates  provided  by  an  V  x  M  image.  A  point  located  at  pixel  coordinates  (x',y')  in 
the  image  frame  can  be  transformed  into  a  position  vector  using  Equation  2.12  [13]. 

^The  two  quadrant  inverse  tangent  is  again  used  in  Equation  (2.8)  to  define  the  length  of  /' 
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(2.12) 


A 

-  0^ 

Y 

= 

Oy-y 

Z 

f 

2.4.3  The  Camera  with  Intrinsic  Parameters 

If  in  addition  to  the  angle  to  the  prineipal  point  and  the  FOV,  the  aspeet  ratio  and  skew 
factor  are  also  known,  the  intrinsic  parameters  or  calibration  matrix  can  be  used  to  solve 
the  astrometry  of  an  image.  Like  the  frontal  pinhole  camera  model,  this  method  defines  the 
pixel  {x' ,y)  where  the  desired  object  is  located  and  then  transforms  it  to  a  position  vector, 
as  shown  in  Figure  2.5. 


Figure  2.5:  Transformation  from  pixel  coordinates  to  \XYZY 


The  intrinsic  parameters  camera  model  given  in  Equation  2.13  can  be  used  to 
determine  the  real  world  position  of  an  object  in  a  CCD  image.  As  long  as  the  intrinsic 
parameters;  aspect  ratio  {if),  skew  factor  (5^),  and  focal  length  (/')  are  known,  this  model 
can  define  the  pixel  location  of  an  object  in  a  normalized  reference  frame  creating  a  position 
vector  [56]. 
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(2.13) 


y  f  sj'  O,  X 

A  y  =  0  Tif  Oy  Y 

1  0  0  1  z 

Where  A  in  Equation  2.13  is  once  again  a  normalizing  scale  factor  and  the  skew  factor  is 
proportional  to  cot(e),  where  e  is  the  angle  between  the  and  image  axes.  In  most 
systems  e  is  very  nearly  90®  so  5^  will  be  assumed  equal  to  zero  [56].  The  inverse  process 
to  determine  the  real-world  location  of  an  RSO  given  the  focal  plane  pixel  location  is  given 
in  Equation  2.14. 

X  x'  -O^ 

Y  =  (2.14) 

t] 

z\[r 

In  addition  to  the  intrinsic  parameters,  to  use  this  approach  knowledge  of  the  extrinsic 
parameters;  and  do  defining  the  boresight  look  angle  of  the  system  is  required  to  correctly 
transform  the  image  to  the  ECI  coordinate  system  [56].  Once  a  position  vector  in  the  ECI 
reference  frame  is  available,  a  spherical  coordinate  transformation  can  be  performed  to 
yield  the  look  angles  from  that  sensor  to  the  object,  successfully  extracting  the  angles  data 
from  the  image  allowing  for  orbital  determination  or  analysis. 

2.4.4  astrometry.net 

Eang’s  recently  developed  astrometry.net  software  performs  image  calibration  for 
astronomical  images  through  geometric  hashing.  This  tool  identifies  the  background  star 
field  and  returns  the  principal  point,  look  angle,  and  the  required  information  to  calculate 
the  EOV  of  the  image,  allowing  for  the  determination  of  the  virtual  focal  length  [13,  56,  58]. 
This  information  allows  for  the  intrinsic  parameters  rj  and  /')  of  Equation  (2.13)  to  be 
determined,  yielding  the  required  conversion  from  a  real-world  position  vector  to  pixel 
space.  The  extrinsic  parameters  (oro  and  do)  are  also  returned  [58]. 
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The  astrometry.net  software  has  previously  been  used  for  SSA  research  in  the  tool 
GEODETICA  at  AFRL  and  by  Moomey  at  AFIT  to  resolve  the  astrometry  of  an  RSO  in  an 
Electro-Optical  image  [59,  60].  The  software  embeds  the  astrometry  WCS  into  the  returned 
fits  file,  using  the  USNO  B  catalog^,  allowing  the  user  to  simply  read  off  the  right  ascension 
(a)  and  declination  (6)  of  any  pixel,  eliminating  the  need  to  perform  complex  calculations 
to  find  the  a  and  6  of  the  satellite  streaks.  Furthermore,  calibration  frames  (dark,  flat,  and 
bias  frames)  do  not  impact  the  result  because  CCD  calibration  is  more  important  for  light 
curve  data  than  astrometry  [33].  As  discussed  in  Section  3.12  astrometry.net  will  also  be 
used  in  this  work  as  covered  in  Appendix  A  to  provide  the  right  ascension  and  declination 
of  the  identified  RSO  streak. 

2.5  Streak  Identification  and  Isolation  Research 

Although  astronomers  use  integration  techniques  to  remove  streaks  from  their 
data,  streak  isolation  in  astronomical  images  for  SSA  has  only  recently  been  seriously 
investigated.  In  the  last  several  years,  multiple  methods  have  been  developed  by  different 
organizations  to  extract  streak  data  from  astronomical  images  [61].  Wallace,  for  example, 
has  developed  an  approach  using  an  algorithm  for  streak  detection  within  a  known 
astronomical  collection  [62].  Oniga  developed  a  feature-detecting  algorithm  that  both 
detects  and  classifies  streaks  as  either  RSO,  aircraft,  or  other  using  decision  trees  with 
approximately  98%  accuracy  [61].  For  continuous  sequences  of  unresolved  imagery, 
Sease  and  Flewelling  developed  a  tool  to  detect  and  track  RSOs  called  GEODETICA. 
The  General  Electro-Optical  DEtection,  Tracking,  Identification,  and  Characterization 
Application  (GEODETICA)  tool  uses  phase  congruency,  edge  detection,  and  a  Kalman 
filter-based  multi-hypothesis  point  tracking  framework  to  process  optical  data  and  provide 
both  orbital  and  light  curve  data  in  an  automated  fashion  [59].  GEODETICA  uses 

*The  USNO  B  catalog  is  an  ECI J2000  referenced  catalog,  properly  placing  all  objects  in  the  image  as  if 
the  had  been  observed  from  the  center  of  the  earth  at  1200  on  1  January  2000. 
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astrometry.net  to  resolve  the  astrometry  of  the  RSO  in  the  image,  returning  the  WCS  right 
aseension  and  deelination. 

Sease  has  demonstrated  that  the  matehed  filter  proeess  ean  be  applied  to  teleseopes 
traeking  near  Earth  objeets  that  appear  to  eurve  in  the  image  plane^  by  first  removing 
eurvature  from  star  streaks,  leading  to  homogenous  streak  geometries  despite  the  loeation 
of  the  star  in  the  image  plane.  The  Sease  proeess  depends  on  re-sampling  the  image  about 
the  axis  of  rotation,  and  returning  the  image  in  either  polar  or  spherieal  eoordinates  to 
provide  an  image  where  all  the  star  streaks  are  uniform  overeoming  the  deerease  in  Signal- 
to-noise  Ratio  (SNR)  eaused  by  the  angular  rate  of  the  sensor  [63,  64]. 

The  most  direetly  applieable  researeh  to  the  presented  work  herein  has  been  published 
by  Levesque,  who  uses  a  series  of  baekground  estimations  and  removal  teehniques, 
followed  by  an  iterative  matehed  filter  to  deteet  RSO  streaks  eompletely  eontained  in  the 
image  frame.  This  method  ean  deteet  Vmag  =  15  using  a  0.35-m  teleseope.  Unfortunately, 
the  matehed  filter  requires  the  use  of  a  TEE,  limiting  this  method  to  eatalog  maintenanee 
and  traeking  of  known  satellites  [65]. 

Eurthermore,  Nikolaev  et  al.  showed  that  eomplete  applieation  of  Levesque’s  model  is 
not  always  required.  Applying  only  a  subset  of  Levesque’s  method  allows  for  the  extraetion 
of  brighter  objeets  from  Eleotro-Optieal  data  [66].  Similarly,  Appendix  A  eovers  how  parts 
of  Levesque’s  method  is  applied  for  streak  extraetion  in  this  work. 

2.6  On-going  Parallel  Research 

2.6.1  Air  Force  Research  Laboratory  Research  Efforts 

Currently  AERL  is  developing  new  MHT  orbital  determination  tools  ineluding  SADIE 
and  Ananke  in  parallel  to  AESPC’s  aequisition  of  the  new  Joint  Spaee  Operations  Center 
Mission  System  (JMS)  [67,  68].  These  efforts  are  striving  to  make  the  best  possible  orbital 

®The  curved  appearance  is  caused  by  the  Earths  rotation,  as  the  Earth  rotates  the  celestial  sphere  remains 
stationary  causing  curved  star  streaks  in  images  with  sufficient  exposure,  or  when  tracking  motion  through 
image  stacks  using  GEODETICA. 
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determination  tool  set  available  to  future  space  operators  and  intend  to  address  many 
of  the  concerns  presented  in  [6].  This  work  is  intended  to  complement  those  ongoing 
efforts  with  a  focus  on  observations  made  by  non-AFSPC  (external)  sensors.  An  orbit 
determination  model,  developed  by  Goff  and  covered  in  Appendix  C,  will  be  used  to  show 
the  applicability  of  the  methods  discussed.  This  model  will  be  limited  to  the  processing 
power  routinely  available  at  AFIT  and  will  not  be  as  robust  as  the  expected  results  of  future 
systems.  The  results  of  the  research  effort  presented  in  Chapter  4  are  expected  to  improve 
when  SADIE  or  IMS  are  used  for  orbital  determination. 

2.6.2  Search  and  Determine  Integrated  Environment 

The  Search  and  Determine  Integrated  Environment  (SADIE)  suite  consists  of 
three  development  efforts,  Auto-Search  and  Determine  (SAD),  covariance  based  track 
association,  and  SAD  like  optical  processing,  designed  to  resolve  UCTs.  Auto-SAD  was 
developed  by  the  Naval  Research  Eaboratory  (NRE)  to  link  together  SAD  and  several 
legacy  tools  requiring  manual  implementation  to  resolve  UCT  [7].  SAD  was  originally 
developed  for  the  Air  Eorce  (AE)  Space  Surveillance  System  (SSS)  fence  processing 
and  requires  position  observations  and  executes  in  a  parallel  computing  environment  to 
determine  orbital  elements  from  a  set  of  unknown  objects  using  a  static  algorithm  and 
fixed  gates  for  correlation  [69].  Auto-SAD  is  an  automated  system  designed  to  replace  this 
manually-intensive  process  [7]. 

2.6.3  Resolving  UCTs  with  Search  and  Determine  Integrated  Environment 

Tracks  remaining  after  the  Special  Perturbations  (SP)  Satellite  Identification  (SID) 

correlation  have  been  performed  are  input  as  UCTs  into  SADIE.  Once  SADIE  has  received 
the  UCT  it  compares  it  against  select  candidates  from  the  space  catalog.  If  that  fails  to 
resolve  the  UCT,  SADIE  will  generate  new  candidate  ranges  from  optical  processing  and 
an  adjustable  statistical  correlation  metric  over  as  many  temporally-separated  observations 
as  possible.  Orbit  determination  is  performed  for  each  hypothesis  and  if  not  all  observations 
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are  used  or  fit  intervals  are  inconsistent  with  sensor  performance  that  hypothesis  is  pruned. 
This  process  is  iterative  as  shown  in  Figure  2.6  [7]. 


SADIE  High  Level  Architecture  and  Flow 
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Figure  2.6:  SADIE  high  level  architecture  [7] 


SADIE  was  recently  tested  with  real-world  data  at  a  DoD  facility  where  it 
demonstrated  98%  correlation  of  objects  with  six  or  more  observations  and  97% 
correlation  with  four  or  more  observation  tracks.  This  simulation  demonstrated  that  SADIE 
could  match  the  outcome  of  the  manually  intensive  legacy  process  [67]. 

2.7  Orbit  Determination 

The  focus  of  the  work  of  Chapters  3  through  4  is  on  the  extraction  of  angles  data 
from  satellite  streaks.  The  presence  of  both  angle  data  and  accurate  time  stamps,  provide 
the  required  data  for  orbital  determination  using  a  number  of  well  understood  methods 
including:  Eaplace,  Gauss,  and  the  Double-R  method  [8,  14]. 
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Orbit  determination  routines  leverage  the  statistics  of  RSO  orbit  propagation  and 
sensor  observation  to  solve  the  least-squares  problem.  These  routines  propagate  the  orbital 
solution  state,  as  well  as  an  estimate  of  the  6-Dimensional  Gaussian  uncertainty  in  the  state 
estimate  covariance.  The  Batch  Least  Squares  filter  and  Extended  Kalman  Filter  (EKE) 
are  commonly  used  estimation  routines  [9,  70].  With  the  increase  of  computational  speed, 
many  favor  the  UKF  as  it  does  not  require  linearization  of  the  dynamics  [71,  72].  The 
UKF  implements  the  Unscented  Transform  (UT)  that  performs  nonlinear  transformations 
using  “sigma  points”.  This  approach  precludes  the  need  for  a  state  transition  matrix  during 
covariance  propagation,  which  is  often  complex  for  higher-order  perturbation  modeling 
and  experiences  relatively  substantial  drift  after  short  time  frames  [8].  Furthermore,  recent 
works  have  shown  that  the  UKF  outperforms  the  EKE  in  accurately  capturing  the  estimated 
error  covariance  [73].  Finally,  since  the  UKF  is  a  numeric  sampling  routine,  one  can 
implement  industry  standard  high-fidelity  orbit  propagators  and  sensor  models  to  estimate 
states  and  covariances.  When  performing  orbit  determination  to  validate  the  usefulness  of 
a  result  in  the  following  chapters,  the  method  provided  in  Appendix  C  will  be  used. 

2.8  External  Sensor  Data 

Astronomical  observations  routinely  collect  serendipitous  satellite  streaks  in  the 
course  of  observing  deep  space  [61,  62,  65].  Typically  these  streaks  are  simply  discarded 
and  without  the  proper  application  of  a-priori  data,  such  as  a  North  American  Aerospace 
Defense  Command  (NORAD)  TEE,  many  of  these  streaks  would  otherwise  simply  become 
UCTs. 

In  the  past  decade  AFRE  has  published  a  considerable  amount  of  research  on  using 
small  commercial  telescopes,  similar  to  university  astronomy  equipment,  to  accurately 
track  RSOs.  These  telescopes  have  demonstrated  an  accuracy  better  than  one  arc-second 
and  the  ability  to  successfully  track  objects  in  Highly  Elliptical  Orbit  (HEO)  orbits 
[74,  75].  Recently  DARPA’s  Orbital  Outlook  program  began  purchasing  time  on  amateur 
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and  university  astronomy  telescopes  to  actively  produce  SSA  data  [76].  Additionally,  as 
covered  in  Section  2.5  multiple  methods  have  been  developed  in  the  last  several  years  by 
other  organizations  to  extract  streak  data  from  astronomical  images  [61,  65]. 

Outside  of  universities,  other  astronomy  sensors  of  interest  include  large  survey 
telescopes,  such  as  Pan-STARRS  and  Large  Synoptic  Survey  Telescope  (LSST),  along 
with  space-based  astronomy  equipment  operated  by  National  Aeronautics  and  Space 
Administration  (NASA).  In  addition  to  being  designed  to  detect  faint  objects  which  also 
routinely  collects  serendipitous  satellite  observations,  each  of  these  sensors  provides  unique 
advantages  [77]. 

Pan-STARRS  for  example  is  co-located  with  the  Maui  GEODSS  telescope  and 
operated  by  United  States  Air  Force  (USAF)  10%  of  the  time  [77].  The  USAF  operational 
time  means  the  telescope’s  operational  and  calibration  status  are  known.  The  proximity  to 
Maui  means  that  even  when  performing  its  astronomy  mission,  predicted  RSO  observations 
can  be  removed  from  the  operational  mission  of  GEODSS.  If  GEODSS  is  performing  fewer 
routine  maintenance  observations,  then  the  additional  time  can  be  used  for  operations  where 
higher  accuracy  is  required,  such  as  those  missions  discussed  in  Sections  2.3  and  2.2.2. 

In  contrast  to  Pan-STARRS,  the  ESST  scope  provides  geographic  separation  from 
all  GEODSS  locations.  This  separation  minimizes  combined  observation  errors  by 
lengthening  the  separation  in  position  vectors  and  providing  better  geometry  to  determine 
radial  direction  measurements  [8].  In  addition,  the  Concept  of  Operations  (CONOPS)  for 
ESST  calls  for  observing  the  entire  western  half  of  the  GEO  belt  once  every  three  days  and 
making  its  observations  publicly  available  instantly  through  the  world  wide  web,  providing 
predictable  observations  of  GEO  satellites  in  addition  to  serendipitous  collects  of  other 
orbital  regimes  [78]. 

Although  the  wide  area  survey  mission,  such  as  those  conducted  by  ESST  and 
Pan-STARRS  are  well  suited  for  providing  SSA,  all  astronomical  missions  can  provide 
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useful  data.  In  addition  to  astronomy  missions  AFRL,  has  several  cold  staring  telescopes 
that  routinely  collect  partial  and  full  serendipitous  satellite  streaks.  If  a  method  for 
extracting  position  data  from  these  streaks  was  available,  this  asset  would  immediately 
become  viable  for  SSA  [25]. 

2.9  Summary 

As  this  literature  review  demonstrates,  there  is  a  need  for  improved  SSA  within  the 
DoD.  Several  AFRL  efforts,  including  SADIE  as  discussed  in  Section  2.6.2,  are  limited 
by  the  number  of  observations  available.  Other  AFRL  initiatives,  such  as  those  discussed 
in  Section  2.2.2.3,  require  considerable  time  from  the  handful  of  high  accuracy  telescopes 
available.  These  needs  can  only  truly  be  met  by  increasing  the  number  of  data  sources 
available. 

In  the  current  fiscal  environment  the  cost  associated  with  building,  validating,  and 
maintaining  a  system  precludes  the  obvious  solution  of  increasing  the  number  of  SSN 
sensors.  Other  programs  such  as  DARPA’s  Orbital  Outlook  program  offer  alternatives 
to  commissioning  new  sensors  by  validating  existing  sensors  and  buying  time  on  them 
for  SSA  purposes.  Astronomical  telescopes  often  collect  serendipitous  SSA  data  while 
performing  scientific  missions  as  discussed  in  Section  2.8  and  these  observations  could 
provide  additional  benefits  to  SSA. 

Any  of  the  above  sources  has  the  potential  to  benefit  SSA  for  the  DoD.  The  presence 
of  potential  error  and  distortion  sources  must  be  addressed  before  the  data  can  be  integrated 
into  a  global  SSA  concept  of  operations,  as  discussed  in  Chapter  3.  Furthermore,  as 
exposure  duration  for  these  sidereal  tracking  systems  is  often  long  enough  to  prevent 
the  capture  of  complete  streaks;  alternative  methods  for  using  partial  streaks  need  to  be 
developed,  such  methods  are  developed  in  Chapter  4. 
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3.  Analysis  of  Angle  Data  Extraction  Techniques 


The  science  of  astrometry  aims  to  accurately  define  positions  of  objects  in  reference  to 
the  celestial  sphere  [33].  Several  different  approaches  exist  for  obtaining  position  data  from 
an  image:  Ma  et  al.  covers  a  considerable  number  in  Reference  7.  This  chapter  focuses  on 
the  comparison  of  three  methods  of  astrometry  for  angle  data  extraction:  affine  transform 
astrometry,  the  pinhole  camera  model,  and  the  intrinsic  parameters  camera  model,  also 
known  as  the  intrinsic  parameters  calibration  model.  Having  defined  each  data  extraction 
method  considered,  see  Sections  2.4.1  through  2.4.3;  it  is  known  that  each  method  has 
its  own  requirements  separate  from  the  image,  and  these  requirements  are  summarized  in 
Table  3.1  [33,  56].  Both  the  pinhole  camera  model  and  the  affine  transform  method  have 
been  successfully  used  in  astrometry  in  the  past  to  track  RSOs  [13,  33]. 


Table  3.1:  Requirements  for  each  astrometry  method 


Method 

Known  Angles 

Focal  Length 

Other 

Affine  Transform 

>  3  (stars) 

N/A 

N/A 

Pinhole  Camera  Model 

ao,  do  (degrees) 

virtual  (px) 

FOV 

Intrinsic  Parameters  Model 

ao,  do  (degrees) 

virtual  (px) 

Syj  ^£7  FOV 

Analysis  was  performed  to  determine  how  each  method  is  impacted  by  distortions 
common  in  optical  tracking  systems.  The  analysis  seeks  to  determine  the  strengths  and 
weaknesses  of  each  method.  The  detailed  analysis  of  each  method  should  provide  answers 
to  the  research  questions: 

1.  What  is  the  level  of  accuracy  achieved  by  each  calibration  method? 

2.  Which  method  best  accounts  for  image  distortion? 
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3.1  Generating  a  Baseline  Image 

A  baseline  image,  see  Figure  3.1,  was  ereated  using  eomputer  simulation.  This  image 
eonsisting  of  a  known  star  field  and  four  RSO  streaks,  will  represent  what  we  would  expeet 
to  see  from  a  system  that  was  free  of  error,  imaging  satellites  with  perfeetly  understood 
positions.  Eaeh  star  and  satellite  streak  in  the  image  is  ereated  by  light  striking  the  CCD  in 
a  manner  that  ean  be  modeled  as  a  point  spread  funetion  (PSF). 
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Figure  3.1:  Simulated  baseline  image 


Mathematieally  an  image  taken  by  an  optieal  sensor  takes  the  form  of  Equation  3.1 

[33]: 
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s{x,y)  =  k®  o  {x,  y)  +  hq 


(3.1) 


where  k,  the  eomposite  PSF  of  the  system,  is  the  initial  impulse  response  of  the  system  to  a 
point  source  including  any  far  field  distortion  effects,  o{x,y)  is  the  undistorted  image,  hq  is 
random  additive  noise  from  the  CCD  itself,  and  5(v,  y)  is  the  distorted  image  created  by  the 
system.  The  PSF  accounts  for  diffraction  effects  of  the  CCD,  the  optics  in  the  system,  and 
the  atmosphere  between  the  system  and  the  light  source  [79].  In  non-coherent  imaging 
systems,  such  as  astronomical  telescopes,  the  image  forming  process  can  be  described 
as  a  linear  system,  meaning  that  superposition  applies,  i.e.  if  two  objects  are  imaged 
simultaneously,  the  resulting  image  is  equal  to  the  sum  of  the  objects  imaged  separately. 
Also,  the  image  of  a  complex  object  can  be  viewed  as  a  convolution  of  the  true  object  and 
the  PSF  [80]. 

This  property  of  linear  systems  allows  for  an  ideal  image  to  be  distorted  in  several 
ways  to  simulate  real-world  distortions  by  convolving  of  the  mathematical  model  of  the 
distortion  with  the  ideal  image.  To  study  the  impact  noise  has  on  the  image,  the  baseline 
image  will  be  distorted  through  computer  simulation.  Each  distorted  image  will  be 
processed  by  all  three  astrometry  methods  and  the  displacement  of  select  pixels  (stars  and 
satellites)  from  the  baseline  image  will  be  measured.  RMSE  across  all  displaced  pixels  will 
be  taken  and  plotted  for  each  distortion  to  measure  the  impact  of  each  type  of  distortion. 

3.2  Error  Sources  that  Impact  Position 

There  are  several  sources  of  degradation  that  impact  all  CCD  imaging  systems,  but 
not  all  of  these  sources  impact  the  recorded  position  of  an  object.  This  research  focuses  on 
the  sources  that  impact  where  the  light  falls  and  is  recorded  in  the  CCD  plane.  These  error 
sources  include  tracking  errors,  lens  effects,  jitter,  and  the  atmosphere.  Each  source,  how  it 
was  analyzed,  and  the  results  are  provided  below. 
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3.2.1  Tracking  Error 

Although  tracking  error  in  images  of  RSOs  can  refer  to  the  orbital  position  error  in 
the  TLE  for  the  satellite,  in  this  chapter  it  is  a  reference  to  error  in  the  tracking  system. 
Considerable  error  in  the  tracking  system  can  cause  a  desired  object  of  interest  to  be  missed 
all  together.  Small  amounts  of  error  however  can  still  result  in  the  capture  of  the  desired 
object,  but  will  introduce  a  bias  in  the  process. 

To  simulate  tracking  error  the  baseline  image.  Figure  3.1,  will  be  processed  using  the 
look  angle  of  the  star  to  the  left  of  Ox  Vir  as  the  target  of  the  captured  observation.  As  Ox 
Vir  is  actually  the  star  being  tracked  in  the  image,  this  incorrect  angle  simulates  an  operator 
incorrectly  identifying  the  right  ascension  and  declination,  thereby  introducing  a  tracking 
error. 

3.2.2  Modeling  the  Effects  of  a  Lens  with  Aberrations 

The  traditional  pinhole  camera  model  does  not  include  a  lens,  and  the  intrinsic 
parameters  approach  assumes  a  thin  lens  [56].  Real-world  systems  often  include  lenses 
and  these  lenses  suffer  from  aberrations  that  reduce  the  quality  of  the  captured  image  by 
the  system  as  rays  from  a  point  source  will  fail  to  fall  on  a  common  focal  point  due  to 
the  lens  through  which  they  pass  [33].  If  the  aberrations  are  significant  enough  then  ray 
optics  approach  is  often  used  to  model  their  effect;  however,  if  the  system  is  near  diffraction 
limited  then  wave  optics  is  the  preferred  method  [81].  As  astronomical  optical  systems  are 
designed  to  produce  diffraction-limited  images,  wave  optics  modeling  will  be  used  here 
[33]. 

The  six  most  common  types  of  aberrations:  defocus  {Wd),  spherical  (Vko4o),  coma 
(Vki3i),  astigmatism  iW22i),  held  curvature  (IT220),  and  distortion  (IT311)  [81].  Each  of  these 
aberrations  has  a  different  effect  on  light  as  it  passes  through  the  lens.  For  example  spherical 
aberration  cause  rays  from  sources  other  than  the  principal  ray  to  fall  ahead  of  or  behind 
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the  focal  point,  where  coma  and  astigmatism  both  affect  the  images  away  from  the  optical 
axis  [33]. 

The  lens  image  simulation  approach  used  here  is  defined  by  Voelz  in  Reference  81 
and  will  only  be  briefly  summarized.  First  aberrations  in  the  lens  causes  space  variations  in 
the  incoherent  imaging  process  meaning  that  the  convolution  used  in  Equation  3.1  cannot 
be  used  to  model  the  effects  of  the  lens.  Therefore,  the  superposition  integral,  provided  by 
Equation  3.2,  is  required  [81]. 


U  {uq,  Vo) 


r>oo 

kJ —OO  kJ —O 


1„  (u,  v)  I  h  (uq,  Vo;  u,  v)  p  dudv 


where  iio,  vo  define  the  ideal  image,  and  the  PSE  given  by  |  h  (uq,  vo;  u,  v)  P  is  defined  by 


Equation  3.3  [81]. 


•/w/-  -  ^  ^ 
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The  value  of  w^p  represents  a  multiple  of  the  wavelength  of  the  light  collected  (A  = 
0.555  pm:  the  peak  of  the  visible  spectrum)  and  was  varied  between  0  x  A  and  10  x  A 
to  account  for  common  lens  aberration  affect  listed  above. 

To  produce  the  simulated  image  each  point  in  the  image  plane  was  selected,  a 
normalized  PSE  corresponding  to  the  point  was  generated  and  applied  to  that  point.  After 
each  point  in  the  image  was  manipulated,  the  output  image  was  generated,  shown  in  Eigure 
3.2.  Eigure  3.3  also  shows  the  impact  of  this  type  of  distortion  by  comparing  the  centroids 
of  several  stars  from  Eigure  3.2  compared  to  the  centroids  of  the  same  stars  from  Eigure 
3.1.  Eigure  3.3  show  a  largely  radial  distortion  pattern  that  one  would  expect  from  a  system 
including  a  lens  and  vignetting.  The  impact  of  another  large  error  source,  the  lens,  is  also 
apparent  in  Eigure  3.2.  Eor  a  detailed  discussion  on  this  modeling  technique  please  see 
Reference  81. 
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Figure  3.2:  Aberrations  effects:  image  Figure  3.3:  Aberrations  effects:  pixels 


3.2.3  Jitter 

The  next  source  that  could  impact  position  of  where  the  light  falls  on  the  cells  is  jitter 
in  the  optical  system.  Jitter,  sometimes  called  “shake,”  represents  movement  in  the  optical 
system  causing  light  to  fall  incorrectly  on  the  imaging  plane.  Jitter  can  be  introduced  by  the 
telescope  mount  as  it  tracks  objects,  by  the  shaking  of  the  platform  the  scope  is  seated  on, 
or  even  by  wind  gusts  shaking  the  telescope.  In  well-designed  optical  systems  the  impact 
of  jitter  should  be  minimal  [82]. 

Jitter  can  be  modeled  as  Gaussian  noise  and  directly  effects  the  optical  systems 
response  to  the  PSF.  Assuming  the  telescope  modeled  here  represents  a  well-designed 
optical  system,  we  expect  the  amount  of  jitter  to  represent  a  small  fraction  of  the  FOV  [82]. 
Therefore  jitter  will  be  modeled  using  the  modulation  transfer  function,  given  in  Equation 
3.4  and  represented  by  k  in  Equation  3.1  [82].  The  spatial  frequency  for  the  jitter  (//)  will 
be  kept  below  ^  of  the  EOV,  and  the  amplitude  of  the  jitter  (crjuter)  will  be  varied  between 
0  and  Convolution  of  the  baseline  image  with  Gaussian  filter,  using  Equation  3.1  and 
ignoring  the  CCD  noise,  effectively  simulates  the  distortion  expected  from  jitter  as  can  be 
see  in  Eigures  3.4  and  3.5  [82]. 
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Figure  3.4:  System  jitter  effects:  image  Figure  3.5:  System  jitter  effects:  pixels 


Figure  3.5  shows  that  the  centroids  have  been  offset  slightly  from  where  they  appear  in 
the  baseline  image.  This  offset  is  expected  from  the  blurring,  visible  in  Figure  3.4,  causes 
elongation  of  several  stars,  thus  changing  the  position  of  the  centroid  of  the  PSF. 

3.2.4  The  Atmosphere 

Atmospheric  turbulence  can  have  a  drastic  impact  on  light  passing  through  it. 
Common  effects  include  scintillation,  best  known  for  causing  stars  to  twinkle,  as  well 
as  the  refraction  of  light,  sometimes  leading  to  multiple  copies  of  a  star  in  an  image, 
known  as  speckles.  In  short  exposures,  less  than  1/4  second,  these  atmospheric  differences 
lead  to  extremely  variable  PSF  magnitudes  in  the  imaging  system.  As  the  PSF  of  a  long 
exposure  experiences  turbulence  equal  to  the  sum  of  turbulence  experienced  in  several  short 
exposures  across  the  same  time  frame,  the  atmospheric  turbulence  in  longer  exposures  can 
be  treated  as  a  single  random  event.  Long  exposures  are  often  preferred  for  dim  object 
detection  because  they  have  a  better  SNR  [33]. 
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The  simulated  exposure  time  used  for  this  phase  of  the  researeh  is  10  seconds.  This 
exposure  length  was  chosen  to  be  longer  than  4-5  seconds,  considered  a  minimum  exposure 
by  some  astronomers,  yet  short  enough  to  contain  complete  satellite  streaks  [61]. 

Because  these  long  exposures  can  be  simulated  as  a  single  random  exposure  and  the 
atmosphere  behaves  like  Gaussian  noise,  simulated  atmospheric  turbulence  was  created 
using  independent  random  Gaussian  noise  with  zero  mean  and  normalized  variance.  This 
Gaussian  distortion  was  created  and  convolved  with  the  baseline  image  using  Equation  3.1 
and  ignoring  the  CCD  noise.  This  noise  effectively  blurred  the  stars  and  reduced  their 
magnitude  in  the  image. 

The  results  of  this  noise  can  be  seen  in  Figure  3.6.  Notice  that  several  of  the  stars 
are  no  longer  visible  as  atmospheric  interference  has  increased  the  noise  floor.  What  is 
not  apparent  in  Figure  3.6,  but  does  show  up  in  Figure  3.7,  is  that  the  atmosphere  has  also 
changed  where  the  light  is  falling  on  the  CCD.  Figure  3.7  show  that  the  centroids  of  the 
stars  have  moved  from  their  true  positions  potentially  impacting  the  perceived  position  of 
objects  the  system  is  trying  to  identify. 


Figure  3.6:  Atmospheric  effects:  image 
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Figure  3.7:  Atmospheric  effects:  pixels 
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3.2.5  Multiple  Error  Sources 

To  get  an  idea  of  how  multiple  noise  sourees  effeet  the  ideal  image  to  aid  in 
determining  the  best  overall  method,  two  noise  sourees  where  applied  to  the  ideal  image 
simultaneously.  The  noise  sourees  where  atmospherie  turbulenee  and  the  effeets  of  lens 
aberrations.  These  noise  sourees  were  ehosen  beeause  they  will  be  present  for  any 
terrestrial  based  eleetro-optieal  system.  The  results  of  this  staeking  of  noise  sourees  is 
displayed  in  Figures  3.8  and  3.9. 
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Figure  3.9:  Multiple  noise  sourees:  pix¬ 
els 


These  figures  show  that  the  impaet  of  eombined  noise  sourees  staek  just  as  Equation 
3.1  indieates  they  would.  Figure  3.9  shows  both  the  loss  of  stars  from  the  atmosphere  as 
well  as  the  shift  in  position  eaused  by  both  the  atmosphere  and  the  lens  aberrations. 

3.3  Real-world  Images 

As  a  final  test,  a  series  of  real-world  images  taken  by  the  TeleTrak  network  loeated 
at  AFIT  were  used  to  determine  how  the  astrometry  teehniques  would  funetion  on  non- 
ealibrated  imagery  from  a  real  system.  The  images  used  eontain  the  major  sourees  of 
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position  error  discussed  above.  They  capture  Intelsat  IR,  as  shown  in  Figure  3.10,  as  it 
passed  through  the  FOV  of  the  scope.  The  scope  is  in  sidereal  tracking  mode,  atmospheric 
and  jitter  effects  can  be  observed  in  the  stars,  and  a  lens  was  present  on  the  system  that  took 
the  image.  The  RSO  also  passed  through  the  system’s  FOV  allowing  for  the  effects  of  the 
radial  distortion  near  the  edge  of  the  image  to  be  captured  in  results  as  well  as  data  from 
the  less  distorted  region  near  the  principal  point. 


Figure  3.10:  A  complete  RSO  streak  captured  by  TeleTrak  (contrast  enhanced) 


3.4  Astrometry  Methods  Performance 

The  overall  performance  of  each  astrometry  method  can  be  seen  in  Figure  3.11  where 
the  RMSE  for  each  method  for  each  noise  source  is  provided.  The  RMSE  reported  in 
Eigure  3.1 1  is  derived  from  the  results  across  all  runs  of  that  method.  This  RMSE  analysis 
allows  for  an  objective  comparison  of  the  performance  of  each  method.  The  RMSE  angle 
and  track  error  is  also  provided  in  Table  3.2.  It  is  evident  that  performing  an  intrinsic 
parameters  calibration  results  in  the  best  position  vectors  when  using  distorted  images. 
Each  individual  analysis  and  the  method’s  performance  is  further  analyzed  below. 
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Figure  3.11:  RMSE  for  each  method  and  noise  source 


Figures  3.12  through  3.17  below  provide  the  results  from  the  five  runs  of  four  RSOs 
for  each  method  to  demonstrate  that  the  results  are  both  consistent  and  representative  of  the 
methods  considered. 

3.5  Baseline  Image 

Because  no  estimation  method  is  perfect  it  is  expected  that  each  of  the  methods  used 
here  has  the  potential  to  introduce  error  into  the  position  of  the  observed  objects  in  the 
frame.  The  baseline  image  was  processed  by  all  three  methods  to  observe  any  inherent  bias 
that  the  methods  introduce.  The  results  can  be  seen  in  Figure  3.12. 

The  position  that  streaks  were  actually  placed  on  the  image,  using  the  process  provided 
in  Section  4.5,  are  also  included  as  “Truth”  for  comparison  purposes.  As  Figure  3.12  shows 
both  the  affine  transform  and  intrinsic  parameters  calibration  astrometry  methods  introduce 
little  to  no  bias  in  the  results.  The  pinhole  model  does  introduce  bias  as  shown  in  Figure 
3.12. 

The  distortion  present  in  the  pinhole  camera  model  is  the  focal  length  distortion  as 
expected  given  the  model.  Inspection  of  Equation  2.11  shows  that  only  horizontal  focal 
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Table  3.2:  RMSE  in  degrees  and  kilometers 


Method 

Error  Source 

Angle  (degrees) 

Track  Error  (km) 

Affine  transform 

Tracking  Error 

0.0009 

0.6623 

Eens  Aberration 

0.0804 

59.1833 

Jitter 

0.0847 

62.3632 

Atmospheric 

0.0871 

64.0734 

Combined  Error 

0.2204 

162.167 

Camera/intrinsic  parameters 

Tracking  Error 

0.0014 

1.0354 

Eens  Aberration 

0.0194 

14.303 

Jitter 

0.0183 

13.5111 

Atmospheric 

0.0048 

3.56397 

Combined  Error 

0.0025 

1.87507 

Pinhole  Camera  Model 

Tracking  Error 

0.1063 

78.1902 

Eens  Aberration 

0.0604 

44.4374 

Jitter 

0.0238 

17.5453 

Atmospheric 

0.0624 

45.9114 

Combined  Error 

0.0513 

37.7804 

length  is  considered  in  the  model.  This  result  is  valid  given  the  model  makes  a  square 
pixel  assumption,  unfortunately  square  pixels  are  not  guaranteed  in  real-world  systems  and 
are  not  present  here.  As  the  pinhole  camera  model  does  not  allow  for  scaling,  some  focal 
length  distortion  is  expected.  This  focal  length  distortion  will  also  be  present  in  other 
results  although  they  may  be  small  by  comparison  to  other  error  sources. 
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Figure  3.12:  RSO  position  offset  introdueed  by  method 


3.6  Tracking  Error  Impact 

The  impaet  of  the  indueed  traeking  error  ean  be  seen  in  Figure  3.13.  Here  a  small 
offset  in  the  true  look  angle  of  the  teleseope  vs.  the  intended  look  angle  leads  to  error  in 
both  right  aseension  and  deelination,  slightly  skewing  the  results  provided  by  two  of  the 
three  methods,  and  drastieally  altering  the  results  returned  by  the  pinhole  eamera  method. 

One  method,  the  affine  transform,  is  virtually  unaffeeted  by  the  traeking  error.  This 
result  is  expeeted  beeause  multiple  stars  are  used  to  determine  the  loeation  of  the  RSOs  thus 
a  small  mistake  in  the  teleseope  pointing  angle  is  eompensated  for.  The  intrinsie  parameters 
ealibration  is  also  able  to  aeeount  for  the  traeking  error,  although  not  as  aeeurately,  as  shown 
in  Figure  3.13. 

The  pinhole  eamera  model  is  the  most  degraded  model  returning  positions  that  are 
off  by  an  RMSE  of  0.11°.  That  0.11°  translates  to  36  are-seeonds,  or  approximately  81 
kilometers  of  along-traek  error,  for  a  geosynehronous  satellite.  The  intrinsie  parameters 
ealibration  returned  results  that  were  offset  by  an  RMSE  of  0.0014°,  translating  to  5.04 
are-seeonds,  or  approximately  one  kilometer  of  traek  error.  Einally  the  affine  transform 
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Figure  3.13:  RSO  position  offset  caused  by  tracking  error 


performs  the  best  with  an  RMSE  of  0.0009°,  approximately  3.24  arc-seconds  of  error,  or 
66  meters  of  along-track  error. 

3.7  Lens  Aberration  Impact  Results 

The  impact  of  lens  aberration  is  shown  in  Figure  3.14.  Here  the  intrinsic  parameters 
calibration  camera  model  appears  to  deal  with  the  lens  aberrations  the  best  with  its  results 
appearing  closest  to  the  true  sources  of  the  light,  offset  by  an  RMSE  of  0.02°.  This  minimal 
impact  is  a  result  of  the  fact  that  the  intrinsic  parameters  account  for  the  displacement  of 
the  PSF  caused  by  the  lens. 

With  a  lens  in  place  both  the  affine  transform  and  the  pinhole  camera  model  methods 
move  away  from  truth.  This  result  is  not  surprising  as  the  affine  method  struggles  with  the 
streak  near  the  edge  of  the  image,  and  the  lens  serves  to  increase  edge  distortion.  Given  that 
the  pinhole  camera  model  assumes  no  lens  it  is  expected  that  the  presence  of  a  non-ideal 
lens  will  impact  the  results  provided  by  the  model  [56].  The  RMSE  analysis  reveals  that 
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Figure  3.14:  RSO  position  offset  introduced  by  aberrations 


the  pinhole  camera  model  outperforms  the  affine  transform  method  with  an  RMSE  of  0.06° 
vs.  0.08°. 

3.8  Jitter  Impact  Results 

Jitter  impact  on  the  output  of  an  optical  system  can  be  seen  in  Figure  3.15.  The  RMSE 
analysis  reveals  that  the  intrinsic  parameters  calibration,  with  an  RMSE  of  0.018°  slightly 
outperforms  the  other  methods.  This  accuracy  is  also  apparent  in  Eigure  3. 15  as  the  intrinsic 
parameters  calibration  results  are  closest  to  truth.  The  affine  transform  method  performs 
the  worst  in  the  case  of  jitter,  with  an  RMSE  of  0.0847°,  and  the  pinhole  results  are  close 
to  the  intrinsic  parameters  calibration  results  with  an  RMSE  of  0.023°. 

3.9  Atmospheric  Turbulence  Results 

The  impact  of  atmospheric  turbulence  is  displayed  in  Eigure  3.16.  Eigure  3.16 
shows  that  an  intrinsic  parameters  calibration  accounts  for  the  effects  of  the  atmospheric 
turbulence  the  best  with  an  RMSE  of  0.0048°.  Pinhole  camera  model  performs  worse  than 
the  intrinsic  parameters  calibration  model  with  an  RMSE  of  0.062°,  and  the  affine  transform 
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Figure  3.15:  System  jitter  RSO  position  offset 


method  handles  the  atmospheric  turbulence  the  most  poorly  of  the  methods  considered  here 
with  an  RMSE  of  0.087°.  The  poor  results  provided  by  both  the  pinhole  camera  model  and 
the  affine  transform  method  is  a  result  of  the  atmospheric  turbulence  shifting  where  the 
light  falls  on  the  CCD  as  shown  in  Figure  3.7. 
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Figure  3.16:  Atmospheric  turbulence  RSO  position  offset 
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3.10  Multiple  Error  Sources 

As  expected,  the  results  of  the  multiple  error  sources  simulation,  shown  in  Figure  3.17, 
bears  similarities  to  the  results  of  the  two  error  sources  used  to  create  it.  The  dominant 
effect  is  similar  to  the  impact  of  the  lens  on  the  image,  but  the  angles  are  spread  out  from 
the  atmospheric  turbulence  experienced. 
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Figure  3.17:  RSO  position  offset  introduced  by  multiple  error  sources 


As  shown  in  Figure  3.17,  the  camera  calibrated  with  intrinsic  parameters  model 
remains  closest  to  truth,  with  an  RMSE  of  0.0025°.  The  pinhole  camera  model  is  the 
next  best,  with  an  RMSE  of  0.05°,  and  the  affine  transform  method  providing  an  RMSE  of 
0.22°. 

3.11  Real-World  System  Images 

The  simulated  results  provide  an  understanding  of  how  each  type  of  distortion  effects 
the  observed  position  of  an  RSO.  Indicating  that  when  the  intrinsic  parameters  are  available 
they  can  provide  the  most  accurate  position.  Using  this  knowledge  to  interpret  the  real- 
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world  data  allows  for  the  determination  of  how  using  another  method  may  affect  the  RSOs 
recorded  position. 

Figure  3.18  shows  the  results  of  performing  the  astrometry  methods  on  a  series  of 
images  taken  by  a  real-world  system.  The  pinhole  and  intrinsic  parameters  camera  models 
return  very  similar  results,  particularly  near  the  images  principal  point,  indicating  that  a 
well  calibrated  system  is  responsible  for  capturing  the  image  [56].  Deviation  is  noticeable 
between  the  pinhole  and  intrinsic  parameters  methods  near  the  edge  of  the  image.  Overall 
an  RMSE  difference  from  the  intrinsic  parameters  method  of  0.00003°  right  ascension  and 
0.0397°  declination. 

The  affine  transfer  method  will  offset  from  these  other  two  methods,  providing  an 
RMSE  difference  from  the  intrinsic  parameters  method  of  0.64°  right  ascension  and  0.8052° 
declination.  The  distortion  in  right  ascension  is  caused  by  a  lack  of  a  good  candidate  star 
on  the  right  side  of  the  image  and  the  noticeable  change  in  declination  is  caused  by  the 
weighted  centroids  of  the  chosen  stars  changing  pixel  location  from  image  to  image.  Even 
without  these  effects,  the  offset  from  the  intrinsic  parameters  return  is  considerable. 
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Eigure  3.18:  Real-world  system  image  detected  position 
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3.12  Summary 

The  analysis  given  above  indicates  an  uncalibrated  image  including  a  streak  created 
by  an  RSO  can  be  processed  to  extract  angles  data.  All  three  methods  explored  return 
simulated  position  data  to  within  0.2°  of  the  true  position  without  prior  calibration.  These 
results  are  representative  of  potential  serendipitous  SSA  sources,  as  most  astronomical 
data  sources  will  provide  images  that  are  well  defined  by  a  star  catalog.  Other  astrometry 
methods  may  also  provide  equally  accurate  results  and  depending  on  the  source  of  the  noise 
a  given  method  may  be  preferred  over  the  others. 

Sease  et  al.  uses  the  tool  astrometry.net,  to  perform  the  astrometry  process  in 
GEODETICA  in  Reference  59.  Lang  provides  a  complete  description  of  the  astrometry.net 
geometric  calibration  in  Reference  58  and  it  outperforms  the  intrinsic  parameters 
calibration  for  real-world  images,  by  approximately  0.5°,  therefore  astrometry.net  will  be 
used  to  process  the  real-world  data  used  in  Chapter  4  as  outlined  in  Appendix  A. 
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4.  Updating  Track  Data  from  Partial  Satellite  Streaks 


This  chapter  examines  the  feasibility  and  utility  of  performing  positional  updates  for 
a  space  object  catalog  using  metric  data  obtained  from  streaks  gathered  by  astronomical 
telescopes.  The  focus  of  this  work  is  on  processing  data  from  three  possible  categories: 
streaks  that  only  enter,  streaks  that  only  exit,  or  streaks  that  cross  completely  through  the 
astronomical  image.  The  methodology  developed  can  be  applied  to  any  dedicated  SSA 
sensors,  or  astronomical  imagery,  to  extract  data  from  serendipitous  streaks  gathered  while 
observing  other  Resident  Space  Objects  (RSOs)  or  astronomical  objects.  The  focus  of 
Phase  II  work  will  be  to  answer  the  following  research  questions: 

1.  Can  observations  of  tracks  only  partially  in  the  image  be  used  to  improve  SSA? 

2.  Can  observations  of  tracks  passing  fully  through  the  image  be  used  to  improve  SSA? 

4.1  The  Challenge 

Traditional  orbital  determination  uses  bounded  streaks,  along  with  image  time  and 
exposure  duration,  to  determine  an  accurate  observation  time  and  develop  position  vectors. 
A  bounded  streak  will  be  defined  in  this  work  to  be  a  streak  with  complete  timing  data, 
provided  by  the  presence  of  both  endpoints  clearly  in  the  image,  as  shown  in  Figure  3.10. 
Identification  of  the  center  of  a  streak  is  used  to  provide  precise  timing  of  the  RSO  [8,  14]. 

Streaks  lacking  one  or  both  end  points,  as  shown  in  Figures  4. 1  and  4.2,  do  not  provide 
the  necessary  timing  knowledge  to  determine  accurate  position  and  velocity  vectors  for  the 
RSO  using  traditional  methods.  An  alternative  approach  addressing  the  lack  of  timing 
knowledge  is  needed  if  partially  bounded  and  unbounded  streaks  are  to  be  used  to  perform 
SSA  functions. 
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Figure  4.1:  A  partial  RSO  streak  captured  by  an  optical  sensor 


Figure  4.2:  An  RSO  streak  passing  through  a  TeleTrak  image 


To  avoid  confusion,  the  term  partially  bounded  streak  will  be  used  to  describe  an  image 
of  an  RSO  with  only  a  single  end  point  visible,  as  shown  in  Figure  4.1,  while  a  unbounded 
streak  will  refer  to  a  streak  that  passes  completely  through  the  image  as  shown  in  Figure 
4.2. 
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4.2  Methodology 

To  determine  the  usefulness  of  partially  bounded  and  unbounded  streaks  two 
experiments  were  conducted  using  observational  data  collected  on  Intelsat  IR.  The  data 
used  was  collected  by  the  TeleTrak  network  as  described  below.  The  raw  data  of  two  second 
exposures  was  used  to  provide  a  control  data  set.  These  data  were  also  stacked  to  provide 
simulated  longer  exposures  and  the  equivalent  streaks  that  longer  exposures  provide.  In 
this  manner,  the  same  observations  could  be  used  to  generate  both  “truth”  positions,  the 
centroid  of  the  RSO  in  a  single  frame,  several  images  of  partially  bounded  streaks,  and 
one  unbounded  streak  for  each  night  of  observations,  to  serve  as  experimental  groups. 
Each  simulated  image  and  corresponding  streak  was  processed.  Comparisons  between  the 
results  for  the  control  group  and  the  experimental  groups  are  made  to  analyze  performance. 

4.3  Hardware  Setup  and  Data  Simulation 

The  TeleTrak  network  has  a  variety  of  optical  and  electromagnetic  sensors.  Originally 
designed  to  track  low  Earth  orbiting  satellites,  the  network  has  recently  started  tracking 
geosynchronous  satellites.  The  data  provided  for  this  experiment  was  a  series  of  two- 
second  exposures,  each  consisting  of  120  stacked  frames  taken  by  an  astronomical  camera 
attached  to  an  80-mm  telescope  on  a  mount  in  sidereal  tracking  mode. 

The  provided  images  were  then  processed  to  extract  the  centroid  location  of  each  RSO 
and  provide  positional  data  using  the  traditional  method,  as  discussed  in  Section  2.2.1.  The 
positional  offset  of  these  streaks  from  the  predicted  position  provides  a  baseline  accuracy 
to  compare  against  the  partially  bounded  streak  data  when  examining  residuals  and  other 
orbital  data. 

4.3.1  Partially  Bounded  Streak  Simulation 

To  simulate  partially  bounded  streaks,  the  two-second  exposure  images  were  also 
stacked  to  create  images  similar  to  Eigures  4.3.  The  stacking  of  the  two-second  exposures 
created  the  simulation  of  images  with  longer  exposure  times  to  provide  partially  bound 
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streaks.  Streaks  with  only  a  single  end  point  were  created  so  that  the  visible  edge  was 
present  in  multiple  places  in  the  image  as  the  satellite  passed  through  the  FOV.  This  process 
allowed  for  data  extraction  from  both  the  center  and  the  edge  of  the  image  to  simulate 
correctly  the  real-world  captures.  Gaps  are  visible  in  the  streak  because  the  system  uses  an 
electronic  and  not  physical  shutter,  the  processing  algorithm  is  designed  to  ignore  a  certain 
amount  of  space  to  see  a  complete  streak. 


Figure  4.3:  A  simulated  partially  bounded  RSO  streak 


4.3.2  Unbounded  Streak  Simulation 

Next,  the  unbounded  streaks  were  created  by  stacking  all  of  the  two-second  exposure 
images  from  a  single  night’s  pass  to  generate  a  streak  completely  through  the  image,  as 
shown  in  Figure  4.4.  Additional  images  of  the  star  field  were  also  stacked  to  lengthen  the 
simulated  exposure  beyond  what  was  required  to  create  the  unbounded  streak.  The  number 
of  additional  images  was  chosen  to  increase  the  length  of  the  exposure  as  much  as  possible, 
without  capturing  a  second  RSO  in  the  image.  Multiple  RSOs  in  the  image  increases  the 
difficulty  of  correctly  correlating  the  objects  to  their  TLE  and  could  lead  to  cross-tagging 
satellites. 
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Figure  4.4:  A  simulated  unbounded  RSO  streak  passing  through  the  frame 


4.4  Assumptions  and  Limitations 

To  determine  correctly  the  position  and  an  associated  time  of  observation  of  an  RSO, 
several  assumptions  were  made.  First,  the  source  RSO  is  known.  Here,  three  nights  of 
observational  data  was  used  to  correlate  the  RSO  to  a  single  satellite  after  the  observation 
was  made.^'^  Second,  the  FOV  is  small  enough  that  the  orbital  arc  can  be  estimated  as 
a  straight-line  distance,  but  large  enough  so  that  astrometry.net  software  can  successfully 
resolve  the  image’s  astrometry.  Finally,  the  length  of  the  exposure,  and  its  stop  time,  are 
known. 

4.5  Predicted  Position 

The  results  of  this  work  are  also  limited  by  the  accuracy  of  the  image  data  provided, 
the  TLE,  and  the  propagator  used.  To  estimate  the  amount  of  error  expected  in  the 
results,  a  prediction  will  be  accomplished  using  the  TLE  retrieved  prior  to  the  first  night 
of  observations.  Eor  these  experiments,  the  TLE  will  not  be  updated  after  the  first  night, 
to  prevent  a  catalog  update  from  skewing  the  results  of  the  experiment.  The  TLE  of  the 

'‘'the  presence  of  the  streak  can  also  be  predicted  in  advance  for  some  survey  missions  such  as  the  LSST 
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observed  satellite  will  be  propagated  using  NORAD’s  SDP4  propagator,  and  Systems  Tool 
Kit  (STK)’s  High  Precision  Orbit  Propagation  (HPOP)  propagator  [83].  Known  sources 
of  error  will  be  estimated  and  modeled  to  provide  a  more  precise  estimated  position.  This 
prediction  will  then  be  used  to  create  residuals  from  the  recorded  observation  as  discussed 
above. 

Even  prior  to  the  first  observation  it  is  possible  to  estimate  the  amount  of  position  offset 
caused  by  the  TEE.  Accuracy  of  TEE  data  for  geostationary  orbits  and  specifically  Intelsat 
satellites  have  been  studied  in  great  detail  in  the  past.  Chan  and  Navarro  compared  the 
publicly  accessible  TEE  with  INTEESAT  proprietary  orbital  elements  determined  using 
their  own  dual  station  ranging  system.  The  purpose  of  their  comparison  was  to  find  out 
how  reliable  the  TEE  was  for  conjunction  detection  analysis.  After  comparing  data  for  a 
full  year  Chan  found  that  the  orbital  differences  showed  an  average  of  30-60  km  difference 
in  position,  when  the  most  recent  TEE  was  propagated  to  the  observational  epoch  [43]. 
Eurther  studies  have  been  performed  on  the  accuracies  of  GEO  TEEs  in  general.  Eriih  and 
Schildknecht  performed  a  four  year  study  and  after  evaluating  13  GEO  objects  found  an 
average  of  25  km  along-track  and  10  km  cross-track  difference  after  propagating  the  most 
current  TEE  to  the  current  epoch.  They  also  found  that  the  error  grows  the  further  from 
the  TEE  epoch,  reaching  hundreds  of  kilometers  of  track  error  in  just  a  couple  of  days 
[44].  Vallado,  published  a  study  in  2007,  comparing  precise  orbital  ephemeris  data  to  STK 
HPOP  propagator,  and  found  it  more  accurately  predicted  the  position  of  GEO  satellites 
[84]. 

Eor  these  experiments,  Intelsat  IR  was  observed  for  three  consecutive  nights  by  the 
TeleTrak  network  of  sensors  located  at  the  Air  Eorce  Institute  of  Technology.  The  epoch 
time  on  the  TEE  used  was  approximately  a  day  before  the  first  expected  observation.  Based 
on  the  work  provided  in  References  43  and  44,  it  is  expected  that  less  than  60  km  of  along- 
track  error  and  less  than  10  km  of  cross-track  error  will  be  present  due  to  TEE  propagation 
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using  SDP4.  For  propagation  using  HPOP  an  anticipated  error  of  less  than  1  km  of  along- 
track  and  cross-track  error  will  be  present,  based  on  the  work  provided  in  Reference  84. 

In  addition  to  the  position  offset  caused  by  the  TLE  propagation  some  discrepancy  in 
where  the  light  focuses  on  the  CCD  is  expected  and  will  have  an  impact  on  the  position  as 
well.  Although  the  optical  system  used  is  diffraction  limited  no  design  is  perfect,  and  some 
distortion  is  expected.  To  determine  the  predicted  position  of  the  RSO  in  the  image,  the 
intrinsic  parameters  camera  model  will  be  used.  The  method  is  based  on  the  camera  model 
given  in  Equation  (4.1).  The  camera  model  is  designed  to  predict  the  pixel  location  of  an 
object  on  a  focal  plane,  given  its  real-world  position  [56,  85]. 


f  sj'  O, 
0  r^f  Oy 
0  0  1 


R  t 

0  1 


(4.1) 


To  use  this  model  to  predict  camera  performance  requires  knowledge  of  the  system’s 
intrinsic  parameters  (/',  i’e,  t].  Ox,  and  Oy),  provided  in  Table  4.1,  as  well  as  the  extrinsic 
parameters,  to  define  R  and  t  the  required  rotation  matrix  and  translation  vector  to  transform 
the  position  vector  from  real-world  coordinates  to  pixel  space  coordinates,  as  discussed 
below.  Additionally  a  Cartesian  coordinate  position  vector  for  the  real-world  location  will 
need  to  be  provided. 


Table  4.1:  Intrinsic  and  Extrinsic  Parameters  for  the  Optical  System 


Intrinsic  Parameters  (px) 

f’ 

Ox 

0, 

M 

N 

Value 

62888 

0 

450.5 

600.5 

1 

900 

1200 

Most  star  catalogs  provide  right  ascension  and  declination  angles  for  the  positions  of 
stars  and  most  satellites  define  their  position  with  the  most  recent  TEE.  Stellarium  does  an 
excellent  job  of  converting  the  TEE  to  a  right  ascension  and  declination  look  angle  at  the 
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required  time,  using  the  SDP4  propagators  [86].  When  this  method  is  used,  conversion  to 
Cartesian  direction  can  be  accomplished  using  Equation  4.2. 

X  cos(6)  cos(q') 

Y  =  cos(5)  sin(Q:)  (4.2) 

Z  sin(d) 

Alternatively,  STK  can  be  used  to  provide  a  predicted  Cartesian  coordinate  vector,  using  a 
number  of  different  propagators,  as  was  done  for  the  unbounded  streak  experiment.  Once 
these  Cartesian  coordinates  are  provided.  Equation  4. 1  can  be  used  to  predict  pixel  location. 

The  extrinsic  transformation  required  in  Equation  4.1  is  to  assure  alignment  of 
the  Cartesian  coordinate  vector  and  the  camera  model’s  reference  frame.  Using  the 
known  principal  point  of  each  image,  an  Euler  angle  rotation  will  provide  the  required 
transformation,  first  rotate  about  the  .fiT-axis  by  ao  -  90°  to  align  the  A-axis  with  the  I  axis 
and  the  T-axis  with  the  J.  Then,  rotate  about  the  A-axis  by  do  “  90°  to  bring  the  K  and 
Z-axis  into  alignment.  This  rotation  is  accomplished  using  Equation  (4.3)  [8]. 

cos(ao  -  90°)  sin(ao  -  90°)  0  10  0 

R  =  -siniaQ  -  90°)  cos(ao  -  90°)  0  0  cos(6o  -  90°)  sin(6o  -  90°)  (4.3) 

0  0  10  -sin(6o- 90°)  cos(So  -  90°) 

To  verify  the  model,  a  star  field  of  20  stars  with  known  look  angles  was  modeled.  The 
angles  were  processed  to  predict  the  pixel  locations  of  the  known  stars,  given  the  intrinsic 
and  extrinsic  parameters  of  the  camera  used.  The  output  was  verified  by  comparing  the 
pixel  locations  returned  to  the  locations  of  the  stars  in  a  real  image  taken  by  the  camera. 
The  results  of  the  camera  model  verification  are  shown  in  Eigure  4.5,  predicted  position 
overlying  the  image  positions,  for  visual  comparison.  The  overlay  of  the  true  position  on 
the  returned  values  is  almost  an  exact  match,  but  residuals  are  evident  in  the  image.  The 
RMSE  of  the  residuals  provides  a  position  offset  of  2.73  pixels  in  the  x  direction  and  3.4 
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pixels  in  the  y  direetion.  Given  that  the  camera  averages  3.28  arc  -  seconds  per  pixel  this 
offset  translates  to  roughly  two  kilometers  of  orbital  position  error  at  the  GEO  belt. 
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Figure  4.5:  Camera  model  predicted  star  positions  overlying  true  positions  (px) 


The  last  potential  sources  of  position  offset  are  atmospheric  distortion  and  the  time 
data  embedded  in  the  header  of  the  image.  The  impact  of  the  atmosphere  is  well  understood 
and  can  be  modeled  as  white  Gaussian  noise,  with  the  mean  atmospheric  impact  on  an 
optical  sensor  of  approximately  200  m  [44].  The  time  stamp  for  each  image  is  provided 
by  the  USNO  and  is  applied  to  the  image,  by  the  network,  after  integration  completes. 
Discrepancies  in  the  time  stamp  directly  result  in  position  discrepancies.  Recent  calibration 
studies  on  the  equipment  used  have  found  that  the  time  data  can  offset  the  along-track 
position  by  as  much  as  10  arc  -  seconds  [60]. 

Given  the  potential  error  sources  discussed  above,  it  is  expected  that  a  position  offset 
for  the  observation  of  the  target  RSO  will  occur.  Approximately  100  km  along-track  and 
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50  km  cross-track  will  result  as  shown  in  Table  4.2.  The  majority  of  this  error  is  caused  by 
the  propagation  of  the  TLE  and  grows  as  we  move  further  from  the  TLE  epoch  time.  Some 
error  is  expected  from  the  time  meta-data  and  the  remainder,  constant  error,  is  caused  by 
the  accuracy  of  our  camera  model. 


Table  4.2:  Predicted  position  RMSE  for  the  Intelsat  IR  TEE 


SDP4 

HPOP 

Time 

Sensor 

Atmosphere 

SDP4  Total  (km) 

HPOP  Total  (km) 

along-track 

60 

0.8 

2 

1.83 

0.20 

64.03 

4.83 

cross-track 

10 

0.8 

0 

2.28 

0.20 

12.48 

3.28 

4.6  Extracting  Streak  Data 

To  determine  the  real-world  position  of  the  RSO  present  in  the  image,  a  two 
step  process  was  used.  The  streak  pixel  locations  are  extracted  as  described  below. 
Then,  to  determine  the  real-world  location  of  an  RSO  given  the  pixel  location  provided, 
astrometry.net  is  used.  Eor  the  detailed  step  by  step  image  processing  method  used  see 
Appendix  A. 

4.6.1  Streak  Extraction 

Determination  of  the  pixel  location  of  the  streak  often  requires  detection  of  a  streak  at 
or  near  the  background  noise  threshold.  To  achieve  this  low  SNR  detection  a  modification 
of  Eevesque’s  method  was  used.  Primarily,  first  the  background  was  modeled  and  removed, 
next  the  bright  stars  identified  and  removed,  followed  by  a  Hough  transform  performed 
as  described  in  Reference  65.  To  ensure  the  Hough  transform  performed  correctly,  the 
processing  software  was  instructed  to  fill  in  the  dark  lines  caused  by  the  stacking  of  images 
and  treat  the  entire  simulated  streak  as  a  single  line.  Applying  Eevesque’s  full  method 
would  decrease  the  required  visual  magnitude  of  an  RSO  for  detection,  but  was  not  required 
for  the  streaks  in  the  data  used  here. 
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4.6.2  astrometry.net 

To  resolve  the  astrometry  of  the  objects  in  the  image,  an  astronomical  calibration  of  the 
image  was  performed.  The  tool  developed  by  Lang,  commonly  known  as  astrometry.net 
was  utilized  to  perform  this  function  [58].  This  tool  identifies  the  background  star  field 
and  returns  the  principal  point,  look  angle,  and  overlays  the  WCS  onto  the  image.  The 
presence  of  WCS  information  in  the  header  provides  the  pixel  space  to  a  real-world  position 
transformation  to  ECI  coordinates  and  6^^),  allowing  the  pixel  location  of  objects  of 
interest  in  the  image  to  provide  angles  data  required  for  orbital  information  to  be  derived 
[58]. 

4.7  Partial  Streaks 

Partially  bound  streaks  will  be  generated  both  entering  and  exiting  the  image.  Both 
leading  edges,  from  streaks  entering  the  image,  and  trailing  edges,  from  streaks  exiting  the 
image  will  be  considered.  The  direction  of  travel  of  the  RSO  will  be  determined  using 
the  predicted  position  provided  by  the  TLE  as  described  in  Section  4.5.  The  prediction 
provided  will  also  allow  for  the  determination  of  residuals  by  comparing  predicted  to 
observed  positions. 

These  residuals  can  allow  for  the  determination  of  track  error  present  in  the  prediction 
by  identifying  the  difference  between  the  prediction  and  the  observation.  It  is  expected  that 
residuals  caused  by  the  TEE  will  grow  as  the  satellite  gets  further  from  the  Epoch  time 
provided.  Two  consecutive  nights  of  data  will  be  considered  to  determine  what  role  the 
expected  increase  in  the  magnitude  of  the  residuals  will  have  when  considering  partially 
bounded  and  unbounded  streaks  for  the  determination  of  position.  Residuals  caused  by  the 

"The  right  ascension,  a,  is  defined  to  be  measured  positive  to  the  east  in  the  plane  of  the  equator  from  the 
direction  of  the  vernal  equinox;  although  it  can  be  reported  from  0*  -  24^  in  this  work  degrees,  0°  -  360",  are 
used  [8] 

"The  declination,  6,  is  measured  northward  from  the  equator  (0"  -  90")  [8] 
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difference  between  detected  pixel  locations  will  also  be  determined,  as  discussed  below,  so 
that  an  accurate  comparison  of  the  methods  can  be  made. 

4. 7. 1  Determining  Residuals 

To  use  images  in  which  only  single  end  points  of  streaks  are  visible,  the  end  point 
must  be  associated  with  the  start  or  stop  time  of  the  image  exposure.  The  challenge  is 
to  determine  the  exact  location  of  the  end  point:  the  edge  detection  method  suggested  in 
Reference  65  does  not  catch  the  very  beginning  or  end  of  a  streak,  as  shown  in  Figure  4.6. 
Therefore,  the  pixel  identified  as  the  endpoint  of  the  partially  bounded  streak  will  have  an 
incorrect  time  when  associated  with  the  image  time  provided  by  the  meta-data.  This  lack 
of  correlation  results  in  decreased  RSO  position  accuracy.  This  effect  can  occur  at  both 
ends  when  this  method  is  used  on  bounded  streaks. 


The  image  processor's 
detected  edge 


The  actual  edge  corresponding 
to  the  beginning  time  of  the 
exposure 


The  center  of  a  streak  is  associated 
with  the  exposure  center  time 


Figure  4.6:  The  notional  difference  between  a  streak’s  centroid  and  detected  end  point 


To  determine  the  offset  value,  a  position  will  be  derived  using  the  centroid  of  a 
bounded  streak.  A  second  position  will  be  found  using  the  detected  end  point  of  the 
partially  bounded  streak  that  ends  with  the  same  frame  used  for  the  bounded  frame.  Both 
positions  will  be  compared  to  the  predicted  position  provided  by  the  TLE.  The  difference 
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in  the  residuals  represents  the  loss  of  accuracy.  If  this  difference  is  small,  then  the  use 
of  partially  bounded  streaks,  with  a  single  end  point  in  the  image,  is  a  viable  method  for 
updating  position  data  in  a  catalog.  The  RMSE  of  the  residuals  obtained  across  all  of  the 
streaks  will  be  used  to  determine  the  relative  impact  of  only  using  the  partially  bounded 
streak. 

4.8  Unbounded  Streaks 

The  three  nights  worth  of  data  on  Intelsat  IR  were  also  used  to  generate  unbounded 
streaks.  To  determine  the  usefulness  of  the  position  data  in  a  unbounded  streak  orbital 
estimation  was  required,  see  Appendix  C.  The  control  group  was  used  on  its  own  to  create 
orbital  data.  Unbounded  streaks  were  then  substituted  for  each  night  in  turn  to  determine 
the  impact  of  a  unbounded  streak  on  the  orbital  data  for  each  of  the  first  two  nights.  As 
shown  in  Table  4.3,  metrics  for  both  the  extracted  angles  and  the  orbits  were  recorded  for 
analysis. 


Table  4.3:  Simulated  Data  Sets  from  Collected  Imagery  of  Intelsat  IR 


Night  1 

#  Images 

Exposure  (sec) 

Night  2 

#  Images 

Exposure  (sec) 

Night  3 

#  Images 

Exposure  (sec) 

Control  Group 

8 

4-6 

8 

4-6 

8 

4-6 

Data  Set  2 

a 

0 

- 

8 

4-6 

8 

4-6 

b 

8 

4-6 

0 

- 

8 

4-6 

Data  Set  3 

a 

1 

186 

8 

4-6 

8 

4-6 

b 

8 

4-6 

1 

181 

8 

4-6 

The  direction  of  travel  of  the  RSO  (through  the  image)  was  determined  using  the 
predicted  position  from  the  TLE.  The  prediction  also  allowed  us  to  determine  residuals  by 
comparing  predicted  to  observed  positions.  The  full  state  covariance  of  the  satellite  and  the 
Mahalanobis  distance  ('T)  were  tracked.  These  metrics  were  used  to  assess  the  usefulness 
of  the  unbounded  streaks. 
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4.8.1  Determining  Possible  Position  and  Time 

Since  the  RSO  is  known,  the  veloeity  for  a  eireular  orbit  ean  easily  be  determined 
using  Equation  (4.4) [87]: 


V  =  — 


(4.4) 


With  the  veloeity  of  the  satellite  determined,  the  ealeulation  of  the  distanee  traveled  by  the 
RSO  in  pixel  spaee  for  a  given  exposure  length  is  a  multi-step  proeess.  The  first  step  is 
to  ealeulate  the  distanee  traveled  in  kilometers,  using  Equation  (4.5).  Then,  eonvert  that 
distanee  to  radians,  using  Equation  (4.6),  and  finally  arc  -  seconds,  using  Equation  (4.7). 


d  =  vt 


(4.5) 


6  =  tan 


(4.6) 


e 


4.8481368111x10-6- 


rad 


(4.7) 

arc— seconds 

To  eonvert  to  pixel  spaee,  proeessing  the  image  through  astrometry.net  returns  several 
ealibration  variables  ineluding  the  number  of  arc  -  seconds  found  in  eaeh  pixel  whieh  is 
provided  in  the pixscale  field.  With  this  knowledge.  Equation  (4.8)  allows  for  determination 
of  the  distanee  traveled  by  the  RSO  in  pixel  spaee  during  the  time  of  the  exposure: 


dr 


(4.8) 


pixscale 

Onee  the  distanee  traveled  by  the  RSO  has  been  determined,  the  loeation  of  the  motion 
also  needs  to  be  determined.  The  original  simulated  image  (see  Eigure  4.4),  ean  help  bound 
the  problem.  By  definition,  the  unbounded  streak  touehes  two  edges  of  the  image.  This 
knowledge  provides  both  a  No  Earlier  Than  (NET)  and  No  Later  Than  (NET)  time  of 
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observation.  The  NET  point  corresponds  to  the  RSO  starting  the  exposure  just  inside  the 
image  frame  and  moving  the  distance  dpx  in  the  direction  of  travel,  as  shown  in  Figure  4.7. 
Determining  the  time  required  for  the  RSO  to  reach  the  center  of  the  streak  observed  in  the 
image  provides  one  possible  point  in  time  to  bound  the  RSO’s  position. 

The  NET  point  corresponds  to  the  exposure  ending  immediately  after  the  RSO  reaches 
the  opposite  edge  of  the  frame  from  the  NET  point,  having  moved  the  distance  dp^  in  the 
direction  of  travel,  as  shown  in  Figure  4.8.  Determining  the  required  time  for  the  RSO  to 
move  from  the  center  of  the  frame  to  this  edge  provides  the  other  time  bound  needed  to 
solve  the  problem. 

Time  of  the  centroid  point  is  important  as  that  will  be  the  single  observation  used  for 
the  orbit  update  when  unbounded  streaks  are  included.  The  NET  and  NET  time  bounds 
and  the  time  between  them  represent  possible  times  that  the  RSO  may  have  occupied 
the  centroid  location  to  determine  the  best  possible  match  using  a  global  nearest  neighbor 
approach,  as  discussed  below. 


Figure  4.7:  A  notional  NET  streak  Figure  4.8:  A  notional  NET  streak 
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4.9  Results  and  Analysis 

The  results  of  the  experiments  presented  in  Seetions  4.3.1  and  4.8  are  given  below. 
The  data  output  of  each  experiment,  the  analysis  of  that  data,  and  the  determination  of  the 
feasibility  and  utility  of  performing  positional  updates  for  a  space  object  catalog  using  the 
metric  data  obtained  from  the  streaks  gathered  are  provided. 

4.9.1  Partial  Streaks  Results  and  Analysis 

4.9.1. 1  Position  Accuracy  when  using  SDP4 

To  determine  the  standard  deviation  and  the  bias  of  each  observation  error,  an  initial 
calibration  was  performed.  Beginning  with  a  TLE  solution  one  day  prior  to  the  first  night  of 
observations  and  the  set  of  bounded  observations,  the  TLE  is  propagated  to  the  time  of  each 
observation.  Then,  the  states  are  transformed  to  a  predicted  observation  and  the  residuals 
determined  for  both  a  and  6  of  the  centroids  of  the  bounded  group.  After  all  the  residuals 
are  collected,  a  Gaussian  distribution  is  fit  to  all  the  residuals  of  both  observations  using  a 
maximum  likelihood  estimation.  The  mean  of  the  distribution  is  used  as  the  bias  for  the 
observed  data  (image).  This  process  determined  a  bias  of  0.0208°  for  right  ascension  and 
0.0016°  for  declination,  with  standard  deviations  of  0.0296°  and  0.0015°,  respectively.  This 
bias  attempts  to  account  for  the  constant  propagation,  sensor,  atmospheric,  and  time  error 
present  in  the  observation  leaving  only  the  randomly  varying  noise  and  the  actual  position 
differences  between  the  propagation  and  observation  present  in  the  residual. 

4.9. 1.2  Determined  Residuals 

The  residuals  from  the  observations  of  Intelsat  IR  are  presented  in  Eigures  4.9  and 
4.10,  showing  the  difference  in  radians  between  the  position  provided  by  propagating  the 
TLE  and  the  observed  position  for  both  bounded  and  the  partially  bounded  streak  for  the 
same  moment  in  time.  The  difference  between  the  predicted  position  and  the  observed 
position  of  the  center  of  the  bounded  streak  is  approximately  517  firad  in  right  ascension. 
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or  approximately  22  km  of  along-track  error,  and  26  jirad  in  deelination,  or  about  1  km 
eross-traek  error. 
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Figure  4.9:  The  determined  right  ascension  residuals 
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Figure  4.10:  The  determined  declination  residuals 


76 


Figures  4.11  and  4.12  plots  the  dilFerenee  between  the  determined  position  of  eaeh 
streak  when  using  a  deteeted  end  point  versus  a  bounded  eentroid  in  the  final  frame  used  in 
the  staek.  As  shown  in  Figure  4.11  the  majority  of  the  differenee  between  the  methods  is 
expressed  in  right  aseension,  staying  below  200  ^rad,  or  10  km  differenee.  The  deelination 
differenee  is  below  90  jxrad,  or  3.8  km  differenee,  as  shown  in  Figure  4.12. 
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Figure  4.1 1:  The  right  aseension  differenee  between  traditional  and  experimental  methods 


Overall  an  RMSE  analysis  of  the  residuals,  provided  in  Table  4.4,  shows  an  offset 
less  than  3  firad  in  right  aseension  and  12  firad  in  deelination  is  apparent  between  the 
use  of  the  eenter  of  the  streak  and  a  deteeted  end  point.  When  eompared  to  the  predieted 
offset  provided  above;  both  the  right  aseension  and  the  deelination  observed  error  aetually 
perform  better  than  expeeted,  but  are  still  in  line  with  past  studies  [43,  44]. 

4.9. 1.3  End  Point  Localization 

These  results  presented  so  far  eonfirm  that  the  predietion  model  is  working  as  expeeted 
and  that  the  differenee  between  bounded  and  partially  bounded  streaks  is  minimal.  To 
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Figure  4.12:  The  declination  difference  between  traditional  and  experimental  methods 


Table  4.4:  RMSE  for  positional  data 


Type 

Right  Ascension  Qirad) 

Declination  (jirad) 

Position  Difference 

Along-Track  (km) 

Position  Difference 

Cross-Track  (km) 

bounded 

516.52 

25.70 

21.78 

1.08 

partially  bounded 

513.80 

34.35 

21.66 

1.45 

Difference 

2.72 

11.65 

0.12 

0.37 

ensure  that  using  the  endpoint  of  the  streak  over  the  centroid  of  a  bounded  streak  is  not 
unnecessarily  introducing  new  error  requires  further  analysis.  To  ensure  the  algorithm  is 
accurately  localizing  the  end  of  the  streak,  a  comparison  was  made  between  the  detected 
endpoint,  in  pixel  space,  and  the  single  frame  with  the  same  associated  time.  The  residuals 
determined  between  the  two  methods  is  provided  in  Figures  4.13  and  4.14. 


The  end  point  detection  tends  to  identifying  a  pixel  near  the  edge  and  end  of  the 
streak,  and  not  the  exact  end  point.  This  offset  is  caused  by  the  algorithm  ignoring  light 
below  80%  peak  intensity  of  the  streak.  To  better  determine  the  accuracy  of  the  end 
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Figure  4.13:  The  difFerenee  in  x’  pixels  between  edge  and  eentroid  loeation 
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Figure  4.14:  The  difFerenee  in  y’  pixels  between  edge  and  eentroid  loeation 


point  loealization,  a  Gaussian  distribution  is  fit  to  all  the  pixel  residuals  using  a  maximum 
likelihood  estimation.  This  proeess  determined  a  bias  oF  -0.47  px  For  x'  and  1.57  px 
For  y ,  with  standard  deviations  oF  9.39  px  and  6.38  px,  respeetively.  This  bias  is  eaused 


79 


by  the  accuracy  of  the  end  point  detection  algorithm  compared  to  the  centroid  algorithm 
used.  The  weighted  centroid  algorithm  used  on  the  individual  frames  selects  the  peak 
intensity  creating  the  difference  in  selected  pixel.  This  difference  is  independent  from  TLE 
propagation  and  introduces  the  1.43  km  along-track  and  0.6  km  cross-track  offset  observed 
above.  As  Figures  4.13  and  4.14  show,  with  the  exception  of  a  single  outlier  (observation 
8),  all  of  the  pixels  selected  by  the  end  point  detection  algorithm  fall  within  a  three  standard 
deviations,  and  most  fall  within  a  single  standard  deviation  of  the  location  of  the  centroid 
provided  by  the  single  frame  measurements.  Observation  eight  in  Figures  4.13  and  4.14 
occur  when  the  RSO  is  near  a  bright  star,  and  the  proximity  to  the  star  clearly  impacts  the 
algorithm  determination  of  the  endpoint. 

4.9. 1.4  Impact  to  Current  Operations 

Although  telescope  operators  state  that  streaks  happen  all  of  the  time,  exact 
quantification  of  the  frequency  of  streaks  is  difficult  to  produce  [25,  65].  A  quick  survey 
of  operating  GEO  satellites  was  performed  to  estimate  how  often  SSA  telescopes  might 
observe  a  serendipitous  satellite.  Of  the  91  operational  GEO  satellites  in  the  western 
hemisphere,  30  are  within  half  a  degree  of  another,  43  have  less  than  one  degree  of 
separation,  and  60  are  within  three  degrees  of  another  vehicle.  Therefore,  for  a  telescope 
with  a  single  degree  FOV  centered  on  a  GEO  vehicle,  there  is  a  33%  chance  that  the  sensor 
will  serendipitously  capture  a  second  satellite.  If  a  GEODSS  sensor,  with  two  degrees  FOV 
was  used,  it  will  have  a  47%  chance  of  serendipitously  capturing  a  second  satellite.  If  a 
wide-area  sensor,  such  as  Pan-STARRS,  was  used  then  there  would  be  a  66%  chance  of 
capturing  another  satellite  in  the  image  as  a  partial  streak  removing  the  need  to  schedule  an 
additional  observation  [30,  88]. 

4.9.2  Unbounded  Streaks 

To  implement  the  UKF,  a  starting  covariance  is  needed  to  accompany  the  TEE  of 
interest.  Since  the  candidate  RSO  is  in  a  GEO  orbit  and  observed  frequently,  a  high 


80 


accuracy  covariance  as  proposed  by  Horwood  is  used  [73].  Define  the  starting  covariance 
in  equinoctial  elements  such  that  Pq  =  A^A,  where 

A  =  diag(0.05  km,  10“^  10-^  10-^  10■^97  x  10“^  rad)  (4.9) 

The  initial  state  is  pulled  directly  from  the  published  TLE. 

For  the  following  analysis,  the  comparison  of  orbital  data  generated  using  three  subsets 
of  the  data  was  required.  The  first  set  of  data  is  created  strictly  from  the  control  group,  to 
provide  a  baseline  analysis.  The  second  data  set  generates  an  orbit  by  excluding  data  from 
the  night  with  the  simulated  unbounded  streak.  The  second  data  set  provide  estimated 
orbital  positions  assuming  that  the  captured  unbounded  streak  was  thrown  out,  and  the 
orbit  was  determined  using  only  the  two  nights  of  data.  The  final  data  set  incorporates  the 
unbounded  streak  data.  These  data  sets  allow  for  analysis  of  the  usefulness  of  the  position 
data  in  the  unbounded  streak.  The  impact  of  using  the  unbounded  streak  vs.  ignoring  the 
streak  and  using  fewer  observations  is  discussed  below. 

4.9.2. 1  Position  Accuracy 

The  unbounded  streak  experiment  used  an  HPOP  propagator  and  UKF  to  aid  in  orbital 
position  determination,  so  a  second  calibration  was  required  to  accurately  account  for 
these  differences.  The  Kalman  filter  is  the  optimal  estimator  when  the  observation  error 
(noise)  is  Gaussian  and  is  perfectly  captured  by  the  observation  noise  covariance  matrix, 
Rnoise-  To  determine  the  standard  deviation  and  the  bias  of  each  observation  error,  an  initial 
calibration  was  performed.  Beginning  with  a  TFE  solution  one  day  prior  to  the  first  night 
of  observations  and  a  set  of  control  data  telescope  observations,  the  TFE  is  converted  to 
equinoctial  elements  and  then  propagated  to  the  time  of  each  observation.  Then,  the  states 
are  transformed  to  a  predicted  observation  and  the  residuals  determined  for  both  a  and 
6  of  the  centroids  of  the  control  group.  After  all  the  residuals  are  collected,  a  Gaussian 
distribution  is  fit  to  all  the  residuals  of  both  observations  using  a  maximum  likelihood 
estimation.  The  mean  of  the  distribution  is  used  as  the  bias  for  the  observed  data  (image) 
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and  the  standard  deviation  of  the  distribution  is  ineorporated  in  the  Rnoise  niatrix.  This 
proeess  determined  a  bias  of  -0.0263°  for  right  aseension  and  0.0172°  for  deelination,  with 
standard  deviations  of  0.0124°  and  0.0014°,  respeetively. 

4.9.2.2  Best  Streak  Observation 

When  proeessing  unbounded  streak  angle  data,  it  is  neeessary  to  assoeiate  a  time  stamp 
with  the  observation  angles.  This  problem  is  solved  using  a  diserete  optimization  seareh 
routine,  over  the  range  of  possible  times,  that  uses  the  UKF  propagation  step  to  determine 
a  metrie  based  on  the  weighted  residual  for  all  possible  time  steps.  The  weighted  residual 
metrie  is  defined  as 


J  = 

Va 

Vs 

(T a 

1 

0-6 

The  best  estimate  for  the  time  step  is  seleeted  as  the  time  step  assoeiated  with  the 
smallest  value  of  J.  Figures  4.15  through  4.17  demonstrate  this  proeess  for  the  first  night’s 
unbounded  streak.  As  the  figures  show  the  best  possible  time  is  assoeiated  with  the  line 
segment  243  whieh  eorresponds  to  08:27:07  Greenwieh  mean  time  (GMT).  With  the  best 
times  assoeiated  with  the  unbounded  streaks,  an  orbital  estimation  ean  be  performed  and 
an  analysis  of  that  orbit  eondueted. 
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Figure  4.15:  First  pass  through  the  optimizer 
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Figure  4.16:  pass  through  the  optimizer 
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Figure  4.17:  Final  pass  through  the  optimizer 


4.9.2.3  Residuals  Analysis 

The  analysis  of  the  estimated  orbits  begins  with  the  eomparison  of  the  estimated 
residuals.  Figures  4.18  through  4.21  compare  the  a  and  5  residuals  from  the  control  group 
to  the  orbital  estimation  created  while  ignoring  the  night  with  a  simulated  unbounded  streak 
(streak  excluded),  and  the  estimation  using  the  unbounded  streak  (streak  included)  from 
the  first  two  nights  of  observation.  Analysis  of  Figures  4.18  through  4.21  shows  very  little 
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difference  between  these  residuals  with  a  RMSE  difference  of  0.00025°  in  right  ascension, 
and  0.00093°  in  declination. 


The  largest  differences  are  observed  in  declination  on  the  second  night  of  observations 
when  the  unbounded  streak  is  observed  on  the  first  night,  see  Figure  4. 19;  however,  once  the 
second  night  of  observations  are  included  the  presence  of  the  unbounded  streak  improves 
the  solution.  This  improvement  is  evident  in  the  third  nights  residuals,  observations  10 
forward  in  Figure  4.19.  Observation  nine  in  Figure  4.19,  which  is  also  observation  one  in 
Figure  4.21,  is  greater  than  three  sigma  from  our  expected  observation.  This  observation 
represents  the  earliest  frame  with  the  RSO  present  on  the  third  night  and  due  to  both 
increased  atmospheric  turbulence,  and  the  proximity  to  the  edge  of  the  sensor’s  FOV  suffers 
from  increased  distortion.  Because  this  increased  distortion  is  present  in  all  three  data  sets, 
including  the  outlier  does  not  impact  the  results. 
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Figure  4.18:  Estimated  a  residuals  of  2”^  and  3''^  night  (F^  night  unbound  streak) 
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Figure  4.19:  Estimated  6  residuals  of  2”^  and  3''“  night  night  unbound  streak) 
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Figure  4.20:  Estimated  a  residuals  of  3™  night  (2"“  night  unbound  streak) 


4.9. 2. 4  Covariance  Analysis 

Figures  4.22  through  4.25  eompare  the  eovarianee  in  mean  longitude  and  semi-major 
axis  using  observations  from  the  eontrol  group  and  the  simulated  unbounded  streaks,  both 
estimations  ignoring  the  night  with  a  simulated  unbounded  streak  (streak  exeluded),  and 
the  estimation  using  the  unbounded  streak  from  the  first  two  nights  of  observation  (streak 
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Figure  4.21:  Estimated  6  residuals  of  3™  night  (2™  night  unbound  streak) 


included).  These  orbital  elements  were  chosen  because  the  orbital  position  in  physical 
space  is  well-represented  by  these  elements. 
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Figure  4.22:  Mean  longitude  covariance 
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Figure  4.23:  Semi-major  axis  eovarianee 
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Figure  4.24:  Mean  longitude  eovarianee 


In  all  eases,  using  the  unbounded  streak  data  drives  the  eovarianee  eloser  to  the  eontrol 
group.  As  the  measurement  errors  have  been  determined  this  improvement  is  expeeted; 
any  observation,  regardless  of  quality,  should  drive  down  the  eovarianee,  although  the 
improvement  doesn’t  neeessarily  eonverge  to  the  true  position.  Therefore,  the  Mahalanobis 
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Figure  4.25:  Semi-major  axis  eovarianee 


distanee  is  required  to  determine  if  this  deerease  in  the  size  of  the  eovarianee  is  aetually  an 
improvement. 

4.9.2. 5  Mahalanobis  Distance  Analysis 

Figures  4.26  and  4.27  show  the  Mahalanobis  distanee  using  observations  from  the 
eontrol  group  along  with  the  orbital  estimation,  whieh  did  not  inelude  the  night  with  a 
simulated  unbounded  streak  (streak  exeluded)  and  the  estimation  using  the  unbounded 
streak  data  (streak  ineluded)  from  the  first  two  nights  of  observation.  Although  the  data 
does  show  that  the  use  of  the  unbounded  streak  data  ean  be  slightly  less  aeeurate  (with  an 
RMSE  value  of  0.2),  80%  of  the  time  the  Mahalanobis  distanee  eonfirms  that  the  orbital  fit 
using  the  unbounded  streak  data  is  equal  to  or  better  than  the  resulting  fit  when  ignoring 
the  unbounded  streak  data.  For  more  information  on  Mahalanobis  distanees  and  how  it  was 
ealeulated  here  see  Appendix  C 
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Figure  4.26:  Night  1  Mahalanobis  distance 
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Figure  4.27:  Night  2  Mahalanobis  distance 


4.10  Summary 

The  experiments  performed  demonstrate  that  for  the  available  real-world  data  used 
here  partially  bounded  and  unbounded  streaks  can  provide  useful  orbital  position  data  for 
RSOs  in  astronomical  imagery.  End  point  localization,  provided  in  Section  4. 9. 1.3,  allows 


89 


for  partially  bounded  streaks  to  be  used.  These  partially  bound  streaks  provide  a  mean 
degradation  of  3.88%  in  right  ascension  and  259.41%  in  declination,  when  compared  to 
the  control  group. 

The  application  of  a  global  nearest  neighbor  time  estimate,  provided  in  Section  4.9. 2.2, 
for  the  unbound  streaks  allows  for  the  reconstruction  of  the  RSO’s  orbit  that  80%  of  the 
time  is  equal  to  or  better  than  the  orbit  reconstructed  while  ignoring  the  unbound  streak, 
as  shown  in  Sections  4.9. 2. 3  through  4. 9. 2. 5.  The  impact  of  including  unbound  streaks 
provides  a  mean  degradation  of  6.2%  in  right  ascension  and  42.69%  in  declination,  when 
compared  to  excluding  the  unbound  streaks  when  reconstructing  the  orbit. 

In  both  experiments  the  percentage  degradation  of  the  cross-track  position  is  worse 
than  the  along-track  degradation  for  two  reasons.  First,  the  cross-track  position  error  is 
smaller,  exaggerating  any  difference.  Second,  the  edge  detector  used  for  pixel  localization 
picks  a  pixel  on  the  edge  of  the  streak  not  at  its  center. 
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5.  Conclusions  and  Future  Work 


The  work  here  in  sought  to  enable  the  use  of  external  sensors  for  SSA  as  diseussed  in 
Chapter  1,  by  building  upon  past  work  in  extraeting  streaks  from  astronomieal  images  as 
found  in  References  61,  65,  89,  and  90.  By  incorporating  Levesque’s  work  in  Reference 
65,  in  series  with  Lang’s  work  in  Reference  58,  and  orbital  determination  techniques  found 
in  References  8  and  91,  partially  bound  and  unbound  streaks  were  used  to  accurately 
reconstruct  an  orbit  of  a  known  RSO  using  real-world  data.  The  findings  support  the 
operational  utility  of  more  sensors,  and  may  enable  secondary  SSA  missions  from  current 
AFRL  sensors  [5,  25,  92].  This  chapter  will  provide  a  breakdown  of  the  conclusions 
of  the  two  phases  of  research,  the  answers  to  the  proposed  research  questions,  and  a 
recommendation  for  future  work  to  build  upon  the  findings  presented. 

5.1  Phase  I  Conclusions 

The  first  phase,  detailed  in  Chapter  3,  sought  to  find  the  best  data  calibration  technique 
for  extracting  time,  right  ascension,  and  declination  data  from  images.  The  best  calibration 
technique  was  found  by  comparing  several  existing  techniques  to  determine  the  sensitivity 
of  each  process  to  the  distortions  found  in  EO  data.  The  focus  of  Phase  I  work  was  to 
answer  the  following  research  questions: 

1 .  What  is  the  level  of  accuracy  achieved  by  each  calibration  method? 

2.  Which  method  best  accounts  for  image  distortion? 

The  analysis  provided  in  Chapter  3  indicates  an  uncalibrated  image  including  a  streak 
created  by  an  RSO  can  be  processed  to  extract  angles  data.  All  three  astrometry  methods 
explored  return  simulated  position  data  to  within  0.2°  of  the  true  position  without  prior 
calibration.  These  results  are  representative  of  potential  serendipitous  SSA  sources,  as  most 
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astronomical  data  sources  will  provide  images  that  are  well  defined  using  a  star  catalog. 
Other  astrometry  methods  may  also  provide  equally  accurate  results,  and  depending  on  the 
source  of  the  noise,  a  given  method  may  be  preferred  over  the  others. 

The  results  of  this  phase  were  compared  to  Sease  et  al.  using  astrometry.net,  to 
perform  the  astrometry  process  in  GEODETICA  in  Reference  59,  published  in  early 
2015.  Lang  provides  a  complete  description  of  the  astrometry.net  geometric  calibration  in 
Reference  58  and  it  outperforms  the  intrinsic  parameters  calibration  for  real  world  images, 
by  approximately  0.5°,  therefore  astrometry.net  is  the  recommended  tool  for  processing 
real  world  images  that  include  enough  stars  for  a  successful  calibration. 

The  RMSE  analysis  results,  provided  in  Figure  3.11,  shows  that  performing  an 
intrinsic  parameters  calibration  outperforms  the  remaining  methods  by  a  range  of  0.04  - 
0.22°  depending  on  the  noise  source,  clearly  demonstrating  that  it  is  preferred  over  the 
other  two  methods  for  serendipitous  streak  harvesting  if  astrometry.net  cannot  be  used.  The 
requirements  for  knowledge  of  the  intrinsic  parameters  might  limit  the  intrinsic  parameters 
methods  use  to  known  systems,  but  the  real  world  results  indicate  that  even  without  the 
intrinsic  parameters  the  pinhole  camera  model  returns  an  almost  equally  accurate  position 
for  diffraction  limited  systems.  The  findings  of  this  phase  clearly  show  that  accurate  orbital 
reconstruction  can  be  performed  using  an  RSO  streak  in  a  distorted  image,  without  applying 
calibration  frames  (dark,  flats,  and  bias  frames). 

5.2  Phase  II  Conclusions 

The  work  of  Phase  II,  detailed  in  Chapter  4,  examines  the  feasibility  and  utility  of 
performing  in-track,  cross-track,  and  orbital  element  updates  for  an  RSO  catalog  using 
metric  data  obtained  from  streaks  gathered  by  astronomical  telescopes.  This  phase  builds 
on  the  findings  of  Phase  I  by  using  the  astrometry  methods  to  determine  accurate  angle  data 
and  reconstruct  orbital  position.  The  focus  of  this  work  was  on  gathering  and  using  three 
possible  data  categories:  streaks  that  only  enter,  only  exit,  or  cross  completely  through 
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the  astronomical  image.  The  focus  of  Phase  II  work  was  to  answer  the  following  research 
questions: 

1.  Can  observations  of  tracks  only  partially  in  the  image  be  used  to  improve  SSA? 

2.  Can  observations  of  tracks  passing  fully  through  the  image  be  used  to  improve  SSA? 

The  data  used  for  this  phase  was  gathered  using  low  cost  Commercial  Off  the  Shelf 
(COTS)  equipment,  located  in  data  Ohio.  The  nature  of  the  equipment  and  its  location 
provides  real-world  data  that  is  less  than  ideal  and  represents  the  lower  end  of  quality  that 
would  be  expected  from  astronomical  scopes. 

5.2.1  Using  Partially  Bounded  Streaks 

As  shown  in  Section  4.9.1,  high  confidence  exists  that  the  method  used  correctly 
isolates  the  end  of  a  streak  and  selects  a  pixel  associated  with  the  position  of  the  observed 
RSO.  Although  the  pixel  selected  is  offset  by  the  threshold  chosen  to  indicate  the  end  of 
the  streak  and  not  the  actual  pixel  location  that  correlates  to  exposure’s  time  stamp,  it  is 
shown  to  be  within  a  single  standard  deviation. 

These  partially  bound  streaks  provide  a  mean  degradation  of  3.88%  in  right  ascension 
and  259.41%  in  declination  respectively  when  compared  to  the  control  group.  This 
difference  is  reflected  in  the  3-12  fxrad  offset  found  when  using  streaks  with  a  single  end 
point  versus  using  a  streak’s  centroid.  The  impact  of  using  a  single  end  point  corresponds 
to  a  difference  of  less  than  0.4  km  of  orbital  positional  discrepancy  for  objects  in  GEO 
orbits.  These  differences  can  be  directly  attributed  to  differences  in  the  detected  pixels 
used  to  determine  the  angles  used  in  the  orbit  reconstruction.  Further,  the  difference  is 
significantly  smaller  than  the  detected  position  difference  of  more  than  22  km  caused  by 
propagation  inaccuracies. 

These  results  support  the  idea  that  using  a  partially  bounded  streak,  can  be  used  for 
a  satellite  position  update,  as  the  observation  provides  a  more  accurate  location  for  the 
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RSO  than  the  propagated  TLE  alone,  even  with  the  offset  caused  by  the  experimental 
method.  The  decrease  in  track  error  provided  exceeds  the  small  offset  the  new  method 
creates,  allowing  for  enhanced  SSA  from  the  use  of  partial  streaks  when  full  streaks  from 
dedicated  RSO  observations  are  unavailable. 

5.2.2  Using  Unbounded  Streaks 

The  combination  of  the  decreasing  covariance  and  the  equally  good  Mahalanobis 
distance  confirm  that  unbounded  streaks  can  provide  useful  position  and  velocity 
information.  The  covariance  could  be  improved  further  by  using  data  with  sub-second 
time  accuracy,  or  gathered  at  a  location  with  less  atmospheric  turbulence,  which  would 
increase  the  resolution  of  the  data.  The  impact  of  including  unbound  streaks  provides  a 
mean  degradation  of  6.2%  in  right  ascension  and  42.69%  in  declination,  when  compared  to 
excluding  the  unbound  streaks  when  reconstructing  the  orbit.  These  results  support  the  idea 
that  using  a  unbounded  streak  can  provide  position  and  velocity  data,  allowing  enhanced 
SSA  from  the  use  of  unbounded  streaks  in  conjunction  with  dedicated  RSO  observations 
to  perform  an  orbital  update. 

In  addition  to  the  potential  efficiencies  in  tasking  GEODSS  sensors,  discussed  in 
Section  4.9. 1.4,  the  ability  to  use  both  partially  bound  and  unbound  RSO  streaks  may 
enable  existing  sensors  to  contribute  to  the  SSA  mission.  Eor  example  AERE  has  a 
number  of  cold  staring  optical  sensor  that  observe  these  types  of  streaks,  and  it  is  expected 
that  ESST  will  also  observe  multiple  satellites  each  night  it  operates  [25,  78].  Both  of 
these  resources  might  provide  valuable  data,  but  because  this  study  was  limited  in  scope 
additional  work  needs  to  be  performed  to  further  study  the  suitability  of  serendipitous 
satellite  streaks. 

5.3  Future  Work 

The  data  used  for  this  work  was  a  single  satellite  (Intelsat  IR)  observed  by  a  single 
optical  system  over  a  three  night  period  and  processed  by  a  time  consuming  method  that 
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still  requires  significant  manual  input.  Follow-on  work  should  be  focused  on  one  of  three 
areas:  automation,  serendipitous  collects,  and  handling  several  observations  from  multiple 
sensors. 

True  automation  of  the  image  processing  method  outlined  in  Appendix  A  would  allow 
real  time  angle  extraction  from  these  images  enhancing  SSA.  Currently  the  threshold, 
minimum  streak  length,  the  size  of  the  gap  (caused  by  staking)  to  ignore,  and  the  number 
of  peaks  to  include  all  have  to  be  set  manually,  see  Section  A.  1.1.2,  if  these  values  could  be 
automatically  detected  each  image  could  be  processed  in  less  than  a  minute.  Alternatively, 
GEODETICA  provides  exactly  the  type  of  automation  required,  but  it  is  not  yet  able  to 
detect  end  points  of  partially  bounded  streaks,  or  process  unbounded  streaks.  Adding  this 
capability  to  GEODETICA  would  make  it  ideal  for  processing  serendipitous  streaks. 

Although  the  experiments  conducted  in  this  work  indicate  that  astronomical  sources 
can  be  of  use  to  global  SSA,  the  data  was  provided  by  AFIT’s  TeleTrak  system.  The  data 
was  collected  while  observing  young  stellar  objects  over  a  three  night  period,  but  future 
work  should  attempt  to  get  real  world  astronomy  data,  or  AFRL  serendipitous  collects  for 
processing. 

Future  work  should  focus  on  confirming  that  these  methods  also  work  for  other  orbital 
regimes.  Streaks  serendipitously  captured  by  optical  systems  have  a  high  probability  of 
being  Low  Earth  Orbit  RSOs.  The  TEE  for  Eow  Earth  Orbit  (EEO)  objects  are  considered 
less  accurate  than  they  are  for  the  GEO  objects  studied  here.  The  same  experimental 
methods  used  here  can  be  applied  to  bright  EEO  objects  imaged  with  rapid  exposure  times 
once  EEO  data  has  been  made  available. 

Serendipitous  streak  data  from  multiple  locations  provides  the  best  opportunity  to 
leverage  external  sensors  and  maximizes  their  impact  on  SSA  operations.  The  increase  in 
awareness  provided  by  the  external  sensors,  as  well  as  the  benefit  of  geographic  separation 
needs  to  be  balanced  against  data  fusion  challenges.  Exploring  the  correct  methods  for 
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integrating  this  data  in  a  timely  manner  will  be  required  before  external  sensor  data  can  be 
used  operationally. 
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Appendix  A:  Image  Processing  Methods 


Each  step  in  the  proeess  of  extraeting  angles  from  an  RSO  streak  used  in  this 
dissertation  are  well  defined  in  other  sourees,  see  Referenees  [8,  33,  56,  58,  90],  but  are 
combined  here  to  create  a  funetioning  angle  extraetion  tool.  This  appendix  will  cover  how 
the  images  in  the  data  set  where  proeessed,  including  a  brief  overview  of  the  eode  used.  The 
overall  proeess  is  provided  in  Figure  A.  1  with  eaeh  step  eovered  in  the  seetions  below.  The 
images  provided  throughout  this  appendix  are  examples  taken  from  the  data  set  deseribed 
in  Chapter  4. 


pixels 


image 


Astrometry.net 


time  stamp 


a,  6 


Figure  A.l:  Image  Proeessing  Method 
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A.l  Data  Reduction 


To  reduce  the  data  provided  by  TeleTrak  the  images  can  be  converted  to  pixel  locations 
using  the  script  provided  in  Listing  A.l 

Listing  A.l:  Script  used  to  reduce  images  to  pixel  coordinates 


%Used  to  reduce  the  images  in  DataAstrometryPaper  to  the  required  pixel 
%locations . 
close  all;  clc; 

%%  Create  the  directory  of  images  to  process 
%  Point  the  processor  at  a  file  of  png  images 
fileFolder=’ C:\Users\Charlie\Desktop\PartialStreaks\25Septl4 ’ ; 
FileZ=dir(fullfile(fileFolder , fits’))  ; 

%assume  a  sperical  earth  for  now,  this  will  be  more  accurate  if  you  use 

%WGS  data  to  perfect  the  radius  to  the  scope  location. 

r_earth=6378 . 135 ; 

for  ii=2 : length(FileZ) ; 

%Process  each  image  returned  to  extract  the  streak  coordinates 
Image=fitsread(fullfile(fileFolder ,FileZ(ii) .name)) ; 
xy_long( : , : , ii)  =  image_processing_d(  Image  ); 

end 


The  script  given  in  listing  A.l  performs  the  data  reduction  by  calling  two  functions. 
Lines  7  and  8  define  a  directory  of  images  to  be  processed.  Lines  9-15  loops  through  the 
images  in  the  directory  to  reduce  them  to  pixel  locations.  Specifically,  line  15  processes 
each  image  using  the  function  described  in  Section  A.  1.1. 

A.1.1  Detect  Streak  and  Extract  Pixel 

In  Listing  A.l  line  12  passes  each  image  to  the  image 4>rocessing  function  provided  in 
Listing  A. 2.  This  function  uses  several  built  in  MATLAB  functions  to  process  the  image  to 
extract  the  required  pixel  locations  for  use  in  the  camera  model  provided  in  Section  2.4.3. 
Figure  A. 2  is  an  example  image  from  the  data  set  and  will  be  processed  in  this  appendix. 

A.l. 1.1  Detecting  RSO  streaks 

The  steps  taken  in  Listing  A.2  are  a  subset  of  the  steps  provided  by  Levesque. 
Primarily,  first  the  background  was  modeled  and  removed,  next  the  bright  stars  identified 
and  removed,  followed  by  a  Hugh  transform  performed  as  described  in  Reference  65. 
Applying  Levesque’s  method  would  decrease  the  required  visual  magnitude  of  an  RSO 
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Figure  A. 2:  Image  (S  jcy)  Returned  by  the  System 


for  deteetion,  but  was  not  required  for  the  streaks  in  the  data  used  here.  Listing  A. 2  is  a 
funetion  that  requires  an  image  and  returns  the  pixel  locations  of  detected  streaks  in  the 
xy -prime  matrix  and  calls  out  the  pixels  of  the  longest  streak  in  xyJong. 

Listing  A. 2:  Process  an  image  to  identify  RSO  streaks 


function  [  xy_long,  xy_prime  ]  =  image_processing(  Image  ) 

%Function  created  by:  C.  Bellows,  AFIT  2014 

%Process  image  to  remove  background,  stars,  detect  streaks  and  extract 
%pixel  locations  using  a  modification  of  Levesque  et.  al.  method  published 
%in  2007.  (cleaned  up  version!) 

%This  function  takes  the  processed  streak  images  of  streaks  and  processes  them  using  an  edge 
%detector  to  return  an  array  of  end  point  pixels  for  each  streak  in  the 
%image . 

%  Inputs : 

%  Image  the  rotated  image  s_xy 

%  Outputs: 

%  xy_prime  the  xy  pixel  locations  of  the  end  points  of  the  strea 
%%step  1:  preprocess  the  image  to  enable  analysis 

%need  to  convert  to  units  as  matlab  struggles  with  fits  file  manipulation. 

Il=mat2gray(Image) ; 

%need  to  convert  to  a  2D  matrix: 
d=size(Il) ; 

l=zeros(d(l, 1) ,d(l,2)) ; 

J(:, :)=!!(:, :,1); 
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figure(23) ;imshow(J) ;title(’s_x_y  2d  gray  scale’); 

%step  2 :  Model  and  remove  image  background 

%remove  all  foreground  using  morphological  opening.  The  opening  operation 
%removes  objects  that  cannot  completely  contain  the  structuring  element. 
background=imopen(l , strel( ’disk’ , 35)) ; 

%view  the  background  approximation  as  a  surface  to  see  where  illumination 
%varies . 
figure (24) ; 

surf (double(backgroundCl : 8 : end, l:8:end))) ,zlim([0  255])  ; 
setCgca, ’ydir’ , ’reverse’) ; 

%subtract  the  background  approximation  image,  from  the  original  image,  the 
%resulting  image  has  a  uniform  background. 

12=3 -background ; 

figure (2 5);  imshow(J2);  title( ’ s_x_y  background  removed’ ) ; 

%step  3:  adjust  the  contrast  by  1%  at  both  the  high  and  low  ends  of  the 
%spectrum. 

33=imadjust(J2) ; 

J3=imsharpen(33) ; 

figure(26) ;imshow(33) ;titleC’s_x_y  contrast  enhanced’); 

%%step  4:  Remove  stars: 

%image  segmentation 
level  =  graythresh(J3) ; 

BW  =  im2bw(J3,  level); 
image=~BW; 

se  =  strel (’disk’ ,6,6) ; 

image_a  =  imerode(image,se) ; 

image2=J3(: , : , 1) ; 

image3=image2 . *image_a; 

image4=imsharpen(image3 ,  ’radius’,  18); 

figure (2 7);  imshow(image4) ;  title(’s_x_y  star  masked’); 

%%step  5:  change  to  binary  and  pass  to  streak  extractor_2 
s_xy=image4>level- . 2 ; 

figure(28);  imshow(s_xy) ;  title(’s_x_y  streak’); 

%step  6:  perform  the  hugh  transform  and  return  the  x’  y’  pixel  location  of  the 
%streaks. 

[  xy_long,  xy_prime  ]  =  streak_extraction_2(  s_xy.  Image  ); 
end 


In  Listing  A. 3  lines  14  through  21  convert  the  image  to  a  two  dimensional  matrix. 
First,  line  16  converted  the  image  into  a  gray-scale  image.  Second,  line  18  determines  the 
height  and  width  of  the  image.  Then,  line  19  creates  a  matrix  of  the  correct  size,  and  line  20 
fills  that  matrix  with  the  correct  numbers.  Finally,  line  21  displays  the  matrix  as  an  image, 
as  shown  in  Figure  A. 2. 

As  with  Levesque’s  method  the  next  step  is  to  model  and  remove  the  background  of 
the  image.  This  removal  is  to  is  accomplished  from  Steps  23-35.  Line  26  removes  the 
foreground  using  morphological  opening.  The  opening  operation  removing  the  objects 
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that  MATLAB  cannot  completely  eontain  the  strueturing  element.  Lines  29-31  displays  the 
baekground  approximation,  given  in  Figure  A. 3.  Line  34  of  Listing  A.3  then  subtraets  the 
baekground  surfaee  from  Figure  A. 2  and  line  35  displays  the  baekground  free  image,  given 
in  Figure  A. 4. 


Figure  A.3:  Mapped  image  baekground 


To  sharpen  the  edges  of  objeets  in  the  image  lines  37-41  adjust  the  eontrast  of  the 
image  by  1%  at  both  the  high  and  low  ends  of  the  speetrum.  Line  39  adjusts  the  low  end 
of  the  speetrum,  using  imadjust.  Then,  line  40  adjusts  the  high  end  of  the  speetrum  using 
imsharpen.  The  results  are  displayed  using  line  41  as  shown  in  Figure  A. 5. 

Onee  the  eontrast  has  been  enhaneed  removing  the  stars  is  aeeomplished  with  lines  43- 
53.  Line  45  determines  the  threshold  separating  the  foreground  objeets  from  baekground 
noise.  Line  46  then  eonverts  the  image  to  blaek  and  white  using  this  threshold  level.  Line 
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Figure  A.4:  Image  after  mapped  background  removed 


s  contrast  enhanced 

xy 
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Figure  A.  5:  The  contrast  enhanced  image 


48  creates  the  required  strel  disk  and  line  49-51  places  the  disk  over  each  star.  Line  52 
sharpens  the  remaining  image  for  display  on  line  53,  as  shown  in  Figure  A. 6. 
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Figure  A. 6:  The  star  masked  image 


Lines  55-57  converts  the  image  to  binary.  Line  56  performs  the  conversion  using  a 
threshold  slightly  below  the  threshold  automatically  detected  by  MATLAB,  to  ensure  the 
streak  remains  visible.  Line  57  displays  the  new  image,  provided  in  Figure  A.7.  Once  a 
binary  image  of  the  streak  is  available  it  is  passed  to  the  streak  extraction  function  discussed 
below,  on  line  61,  which  returns  the  pixel  locations  of  the  detected  endpoints  of  all  streaks 
as  well  as  marking  the  longest  streak  detected 

A.1.1.2  Determine  Pixel  Locations 

Once  the  binary  image  has  been  passed  to  the  streak  extraction  function,  given  in 
Listing  A. 3,  the  edges  can  be  detected  and  a  Hugh  transform  performed  to  extract  pixel 
locations  for  the  streak.  Listing  A. 3  is  a  function,  called  by  Listing  A. 2  that  requires  a 
binary  image  (s_xy)  and  the  original  fits  image,  and  returns  the  pixel  locations  of  detected 
streaks  in  the  xy^rime  matrix  and  calls  out  the  pixels  of  the  longest  streak  in  xyJong.  The 

'^for  the  experiments  described  in  Chapter  4  these  endpoints  where  used  to  determine  the  centroid  of  the 
streak,  except  when  the  end  points  where  required  for  the  experiment,  as  described  in  Section  4.7 
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Figure  A.7:  The  binary  streak  image 


funetion  provided  in  Listing  A. 3  also  marks  the  streaks  on  the  original  image,  the  longest 
streak  is  eolored  red,  and  any  seeondary  streaks  would  appear  green. 

Listing  A. 3:  Finding  the  pixel  loeation  of  image  streaks 


function  [  xy_long,  xy_prime  ]  =  streak_extractionC  s_xy,  II  ) 
%streak_extraction  returns  pixel  locations  for  streaks  in  an  image 
%Function  created  by:  C.  Bellows,  AFIT  2014 

%  This  function  takes  the  processed  streak  image  and  processes  them  using 
%a  Hugh  transform  to  return  an  array  of  end  point  pixels  for  each  streak  in 
%the  image. 

%  Inputs : 

%  s_xy  the  image  file  created  by  image_processing_c .m 

%  II  the  fits  file  returned  by  astrometry.net 

%  Outputs: 

%  xy_prime  the  xy  pixel  locations  of  the  end  points  of  the  streak 

%%  Satellite  streak  detection 

%% - 

%Extract  edges. 

BW  =  edge (s_xy ,’ canny ;  figure(24);  imshow(BW) ; 

%Display  the  Hough  matrix. 

[H, theta, rho]  =  hough (BW) ; 

figure(25),  imshow(imadjust(mat2gray(H)) , [] , ’XData’ , theta, ’YData’ ,rho, . . . 

’ InitialMagnification’ , ’ fit ’ ) ; 
xlabel(’\theta  (degrees)’),  ylabel(’\rho’) ; 
axis  on,  axis  normal,  hold  on; 
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colormapChot) 

%Find  the  peaks  in  the  Hough  transform  matrix,  H,  using  the  houghpeaks 
%£unction. 

P  =  houghpeaks(H, 50, ’threshold’ ,£loor(0. l*max(H(:)))) ; 

%Superimpose  a  plot  on  the  image  of  the  transform  that  identifies  the  peaks 
xy_prime  =  theta(P( : , 2)) ; 
y  =  rho(PC: , 1)) ; 

plot (xy_pr ime , y , ’ s ’ , ’ color ’ , ’ black ’ ) ; 

%Find  lines  in  the  image  using  the  houghlines  function. 

lines  =  houghlines (BW, theta, rho,P, ’FillGap’ ,20, ’MinLength’ ,500) ; 


%Create  a  plot  that  superimposes  the  lines  on  the  original  image. 
d=size(lines) ; 
temp=zeros(2 , 2 ,d(l , 2)) ; 

figure,  imshow(Il),  hold  on,  title( ’Detected  Streaks’); 

max_len  =  0; 

for  k  =  1 : length(lines) 

xy  =  [lines(k) .pointl ;  lines (k) .point2] ; 
temp(:,:,k)  =  [lines(k) .pointl ;  lines(k) .point2] ; 
plot(xy(: , 1) ,xy(: ,2) , ’LineWidth’ ,2, ’Color’ , ’green’) ; 

%  Determine  the  endpoints  of  the  longest  line  segment 
len  =  norm(lines(k) .pointl  -  lines(k) .point2) ; 
if  (  len  >  max_len) 
max_len  =  len; 
xy_long  =  xy; 
end 


end 

%  highlight  the  longest  line  segment 

plot(xy_long(: , 1) ,xy_long( : , 2) , ’LineWidth’ , 2 , ’Color’ , ’red’ ) ;  hold  off; 


%%  get  the  pixel  values  of  interest 
[ix,iy]  =  find(temp  >  0); 
c  =  [ix  iy] ; 

xy_prime  =  zeros (size (temp) ) ; 
for  kk=l : length(c) 
indices=c(kk, : ) ; 
xpos  =  mod(kk-l,2)+l; 
ypos  =  mod(floor((kk-l)/2) , 2)+l ; 
if  ([xpos  ypos]  ~  [1  1]) 
twosquare  =  zeros(2); 

end 

twosquare (xpos,  ypos)  =  temp(indices(l) , indices(2)) ; 
if  (mod(kk,4)  ~  0) 
j  =  kk/4; 

xy_prime ( : , : , j )  =  twosquare; 

end 

end 

xy_prime=xy_prime( : , : , 1 : j) ; 

% 


0/  0/0/  (VO/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/O/ 

/o  /o/o  7o7o7o7o7o7o7o/o/o7o7o7o7o7o7o7o7o/o7o7o7o7o7o7o7o7o 


end 


To  extract  the  streaks  first  the  edges  in  the  image  are  detected.  This  is  done  using  the 
MATLAB  edge  function.  MATLAB  describes  the  edge  function  as  taking  a  binary  image  as 
its  input,  and  returning  a  binary  image  of  the  same  size,  with  I’s  where  the  function  finds 
edges  and  O’s  elsewhere.  The  canny  method,  finds  edges  by  looking  for  local  maxima  of 
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the  gradient,  was  chosen  because  it  out  performed  other  methods  during  the  verification 
process. 

Lines  18-24  determine  and  display  the  Hough  Transform.  Line  19  uses  the  MATLAB 
hough  function,  which  implements  the  Standard  Hough  Transform.  MATLAB  describes  the 
function  as  being  designed  to  detect  lines  using  the  parametric  representation  of  a  line: 

p  =  xcos(6) +ysm(6)  (A.l) 

Where  the  variable  p  is  the  distance  from  the  origin  to  the  line  along  a  vector  perpendicular 
to  the  line,  and  6  is  the  angle  between  the  x-axis  and  this  vector.  Lines  20-24  then  graph 
the  Hough  Transform  matrix,  as  shown  in  Figure  A. 8. 


-80  -60  -40  -20  0  20  40  60  80 

8  (degrees) 


Figure  A. 8:  The  Hough  Transform  of  the  binary  streak  image 


Once  the  Hough  Transform  matrix  has  been  determined  line  28  uses  the  matlab 
houghpeaks  function  to  find  the  peaks  in  the  Hough  Transform  matrix,  defining  both  the 
number  of  peaks  and  the  threshold  those  peaks  most  cross  to  be  considered.  Lines  29-32 
then  superimposes  a  plot  of  these  peak  locations  onto  the  image  of  the  transform. 
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Line  34  then  uses  the  matlab  houghlines  funetion  to  find  the  lines  in  the  image, 
defining  both  the  size  of  gaps  whieh  may  be  ignored  and  the  minimum  line  length.  Both  of 
these  numbers  are  large  here  because  the  function  was  last  used  to  process  unbound  streak 
data  as  described  in  Chapter  4.  Once  the  lines  in  the  image  have  been  identified  lines  36-53 
of  Listing  A. 3  superimpose  those  lines  back  onto  the  original  fits  file  as  shown  in  Figure 
A.9.  The  longest  line  in  the  image  is  identified  on  line  49  and  returned  by  the  function. 
Finally,  lines  55-72  are  used  to  extract  pixel  locations  of  other  streaks  in  the  image  from 
the  array  they  are  located  in  and  place  them  into  the  matrix  xy  .prime  to  also  be  returned  by 
the  function. 


Detected  Streaks 

2014092Sg4frameg0mm^tellacam^lP26241g947.png 


Figure  A.9:  The  identified  RSO  streak 


A.2  Convert  Pixels  to  Spherical  a  and  6 

Once  pixel  location  for  the  streaks  has  been  determined  those  pixels  need  to  be 
transformed  into  usable  position  vectors.  The  first  step  is  to  average  the  end  point  pixel 
location  for  each  streak  to  determine  the  center  location  of  the  pixel,  see  Section  4.7.1 
for  why  this  is  important.  Then  the  image  and  the  pixel  location  (xy  prime),  can  be 
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given  to  astrometry.net.  As  discussed  in  Chapter  4,  astrometry.net  is  used  to  perform  an 
astronomical  calibration  on  the  original  image,  see  Figure  A. 2.  The  tool  identifies  the 
background  star  field  and  returns  the  astronomical  WCS  meta-data  including  the  principal 
point,  pixel  scale,  and  the  rotation  of  the  image[58].  Each  of  these  pieces  of  information 
is  required  to  find  the  correct  spherical  coordinates  of  the  RSO,  and  create  the  WCS 
information  returned  in  the  header  of  the  new  fits  file.  The  WCS  provides  the  pixel  space  to 
a  real-world  position  transformation  to  ECI  J2000  coordinates,  allowing  the  pixel  location 
of  an  RSO  in  the  image  to  directly  provide  the  angles  data  required  for  orbital  information 
to  be  derived.  Eor  a  detailed  verification  and  validation  of  the  astrometry.net  tool  see 
Reference  58. 

A.3  Right  Ascension,  Declination,  and  Time 

Now,  the  correct  right  ascension  and  declination  have  been  returned.  The  time  stamp 
can  be  pulled  directly  from  the  fits  header  and  the  required  information  for  a  B3  report  used 
by  AESPC  has  been  determined.  Angles  determined  using  this  method  are  used  to  provide 
an  orbital  update  using  additional  methods,  see  Chapter  4and  Appendix  C. 
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Appendix  B:  Verification  and  Validation  of  Angle  Extraction  Techniques 


B.l  Transforming  Camera  Model  Position  Vectors  into  the  SEZ  Reference  Frame 

To  perform  the  verification  and  validation  of  the  models  compared  in  Chapter  3  a 
star  field  was  used  to  ensure  proper  application  of  each  method.  To  ensure  an  accurate 
comparison  the  output  of  each  method  was  transformed  into  a  sensor  reference  frame.  For 
the  affine  transform  method  the  reference  star  locations  were  provided  in  the  traditional 
SEZ  coordinate  system,  returning  other  star  positions  in  the  same  coordinate  system.  Both 
camera  models  had  to  be  transformed  from  the  normalized  \XYZY  vector  to  the  SEZ  sensor 
reference  frame  requires  both  a  rotation  to  line  up  the  coordinate  systems  and  a  translation 
to  co-locate  the  system  origins.  Note  that  the  camera  model  reference  frame  is  not  the 
same  for  the  pinhole  camera  model  as  it  is  for  the  intrinsic  calibration  model.  Therefore 
each  system  will  require  its  own  rotation  matrix  as  given  below  in  Sections  B .  1 . 1  and  B .  1 .2. 

B.1.1  Pinhole  Camera  Model  Position  Vector  to  the  SEZ  Reference  Frame 

To  transform  the  \X  Y  ZY  position  vector  into  the  SEZ  frame  requires  rotating  about 
the  Y  axis  by  (90°  -  cIq),  and  the  intermediate  Z  axis  by  the  azimuth  (fio)  [13].  Eigure  B.l 
provides  an  example  of  how  these  systems  might  be  misaligned.  To  align  these  systems 
first  rotate  around  the  Y  axis  using  Equation  B.l  [8]. 


ROT  I  (90°  -  elo)  = 


I  cos  (90°  -  elo)  0  -  sin  (90°  -  cIq)  \ 

0  1  0 
sin  (90°  -  elf)  0  cos  (90°  -  elo) 

This  rotation  brings  the  zenith  and  z-axis  into  alignment  as  shown  in  Eigure  B.3. 


(B.l) 


Rotating  to  align  the  X  axis  with  the  east  direction  and  the  Y  axis  with  the  south 
direction  can  be  accomplished  with  a  rotation  around  Zenith  equal  to  fio  as  demonstrated  in 
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Figure  B.l:  Offset  between  systems  Figure  B.2:  Aligning  the  Zenith/Z  axis 


Figure  B.3.  Mathematically  this  is  done  using  the  rotation  matrix  defined  by  Equation  B.2 

[8]. 


ROT3(fio)  = 


cos  (J3o) 
-  sin  (J3q) 


sin  (fio) 
cos  (fio) 


0  0 


0 

0 

1 


(B.2) 


With  the  final  rotation  complete  the  axis  of  both  systems  are  aligned  as  shown  in  Figure 
B.4.  Therefore  this  2-3  rotation  provides  a  position  vector  p  in  the  SEZ  coordinate  system 
when  applied  to  the  [XYZY  vector  as  shown  in  Equation  B.3. 


s 

X 

E 

=  ROT3(J3o)ROT2(90‘’ 

-  elo) 

Y 

Z 

Z 

(B.3) 


B.l. 2  Intrinsic  Calibration  Model  Position  Vector  to  the  SEZ  reference  frame 
Eigure  B.5  provides  an  example  of  how  the  intrinsic  calibration  model  and  the  SEZ 
coordinate  systems  might  be  misaligned.  To  assure  alignment  of  the  systems  a  1-3  rotation 
is  required.  Eirst  rotate  around  the  x-axis  by  90°  -  el  as  demonstrated  by  Eigure  B.6. 
Mathematically  this  is  done  using  the  rotation  matrix  defined  by  Equation  B.4  [8].  This 
rotation  brings  the  zenith  and  z-axis  into  alignment  as  shown  in  Eigure  B.6. 
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Figure  B.3:  Aligning  East/X,  South/Y  Figure  B.4:  The  aligned  system 


1  0  0 

ROT\{9Q"  -ek)=  0  cos(90^-e/o)  sin(90^-e/o)  (B-4) 

0  -  sin  (90®  -  e/o)  cos  (90®  -  el^) 


Zenith  Zenith 


Figure  B.5:  Offset  between  two  systems  Figure  B.6:  Aligning  the  Zenith/Z  axis 

Rotating  to  align  the  X  axis  with  the  east  direction  and  the  Y  axis  with  the  south 
direction  can  be  accomplished  with  a  rotation  around  Zenith  equal  to  270  +  /3  as 
demonstrated  in  Figure  B.7.  Mathematically  this  is  done  using  the  rotation  matrix  defined 
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by  Equation  B.2  [8].  The  azimuth  used  here  is  a  referenee  to  the  bore  sight  angle  and  not  a 
traeked  satellite. 

With  the  final  rotation  eomplete  the  axis  of  both  systems  are  aligned  as  shown  in  Figure 
B.8.  Therefore  this  1-3  rotation  provides  a  position  veetor  p  in  the  SEZ  eoordinate  system 
when  applied  to  the  [XYZ]^  veetor  as  shown  in  Equation  B.5. 


s 

X 

E 

=  ROT3{fiQ)  ROT  \  {90^ 

-  elo) 

Y 

Z 

Z 

(B.5) 


Zenith/Z 


Zenith/Z 


East/Y 


South/X 


Figure  B.7:  Rotation  about  3''^  Axis  Figure  B.8:  The  aligned  system 


B.1.3  Spherical  Coordinate  Transformation 

Onee  the  position  veetor  is  defined  in  the  SEZ  eoordinate  system  a  simple  spherieal 
transformation  is  required  to  determine  the  look  angle  (azimuth  and  elevation).  First  the 
azimuth  of  the  traeked  satellite  ean  be  ealeulated  using  Equation  B.6  [93]. 


f  =  aretan 


(B.6) 


Note  that  a  two  argument  are-tangent  is  required  to  assign  the  result  to  eorreet 
quadrant.  Otherwise  the  quadrant  will  need  to  be  assigned  manually  using  the  following 
definitions:  ^  ^  first;  ^  ^  seeond;  5  ^  third;  ^  ^  fourth;  to  ealeulate  the  elevation  use 
Equation  B.7  [93]. 
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(B.7) 


el  =  arcsin  — 

\  IpI 

B.1.4  Translating  the  Position  Vector 

The  origin  of  the  SEZ  eoordinate  frame  is  aetually  loeated  inside  the  sensor  or  CCD 
itself,  eoineiding  the  origin  with  the  Xs,  Xy  system.  Therefore  to  make  the  eamera  model 
position  veetor  eoineide  with  the  SEZ  system  the  origin  must  shift  to  align  the  prineipal 
point  with  its  position  in  real  spaee.  This  translation  is  aeeomplished  by  determining  the 
differenee  between  the  reported  angle  of  the  prineiple  point  and  several  other  stars,  and 
the  knowledge  of  the  objeets  true  position  and  then  adding  the  average  differenee  to  the 
spherieal  eoordinates  returned  by  Equations  B.6  and  B.7  as  shown  in  Equation  B.8. 


^translated 

+ 

^hranslated 

el 

Tel 

(B.8) 


B.2  Verification  and  Validation  of  Astrometry  Tools 

To  ensure  that  the  models  work  as  intended  a  brief  verifieation  and  validation  was 
performed.  The  models  are  validated  as  a  part  of  the  verifieation  proeess,  eaeh  model 
should  return  pixel  speeifie  look  angles  given  the  appropriate  inputs.  To  verify  that  the 
models  work  two  verifieation  steps  are  used.  Eirst,  a  white  box  verifieation  is  aeeomplished. 
Seeond,  the  eorreet  azimuth  and  elevation  for  20  stars  with  known  angles  were  proeessed 
and  the  output  eompared  to  truth. 

B.2.1  White  Box  Verification  and  MATLAB  Implementation 

Implementing  the  mathematieal  theory  outlined  in  Seetions  2.4.1,  2.4.2,  and  2.4.3 
will  be  performed  using  MATEAB  seripts  and  funetion  ealls.  To  ensure  these  seripts  and 
funetions  perform  as  expeeted,  a  white  box  verifieation  of  eaeh  seript  and  funetion  will  be 
performed.  This  verifieation  will  identify  eaeh  equation  in  the  seetions  listed  above  to  a 
line  of  eode  in  the  listings  given  below. 
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The  scripts  for  each  method  provided  in  the  sections  listed  above  will  be  given  in 
Section  B.2.1.1,  B.2.1.2,  and  B.2.1.3.  Each  function  called  by  these  scripts  is  provided  in 
Section  B.2.2. 

B.2.1.1  MATLAB  Affine  Transform  Astrometry  White  Box  Verification 
The  method  for  performing  an  astrometry  using  an  affine  transform  is  provided  in 
Section  2.4.1.  This  method  is  implemented  in  MATLAB  using  the  script  provided  in 
Listing  B.l,  given  below. 

Listing  B.l:  Affine  transform  modei  source  code 


%  This  script  is  designed  to  allow  for  affine  transform  processing  of  an 
%  image: 

%Note:  you  will  need  to  load  a  table  of  the  extracted  pixel  locations  for 
%objects  of  interest. 

%Note:  the  program  will  return  a  look  angle  in  the  same  coordinate 
%system  it  is  provided  reference  angles. 

%%  Given: 

F0V=[1 .402*pi/18Q  1 . ®57*pi/18Q] ;  %degrees  to  radians 
boresight=[214. 1906*pi/18Q  38 . 241944*pi/180] ; 

%beta  stars=213  34’ 51"  213  44’IS2"  214  46’ 37" 

beta_stars=[213 . 580833*pi/18Q:  213 . 733889*pi/180;  2 14 . 776944*pi/180] ; 

%el  stars=38  25’32"  37  38’40"  38  47’37" 

el_stars=[38.425556*pi/188;  37.644444*pi/188;  38 . 793611*pi/188] ; 
ref_pixels=[29. 612983225886452  2 . 364193548387897e+82 ; . . . 

93 . 878129878129878  5 . 983246753246754e+82 ; 5 . 461363636363636e+82  66] ; 

%%  Perform  an  affine  transform  to  produce  the  observed  angles: 
az_el  =  affine_transform(  boresight(l , 1) ,boresight(l , 2) ,  beta_stars, . . . 
el_stars,  ref_pixels,  streak_pixeis) ; 


In  the  script  provided  line  10  is  used  to  provide  the  systems  LOV.  Line  11  provides 
the  bore  sight  azimuth  and  elevation.  Lines  13  and  15  provide  the  azimuth  and  elevation  of 
the  reference  stars  in  the  image.  Lines  16-17  provide  the  pixel  values  for  those  reference 
stars.  The  “affineJransform”  function,  called  on  line  20,  is  provided  in  Section  B.2.2. 1. 

B.2.1.2  MATLAB  Pinhole  Camera  Model  White  Box  Verification 
With  the  understanding  of  how  to  use  the  pinhole  camera  model,  provided  in  Section 
2.4.2,  the  process  is  implemented  MATLAB  using  the  script  provided  in  Listing  B.2. 
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Listing  B.2:  Pinhole  camera  model  source  code 


%  This  script  performs  several  calculations  and  function  calls  to  run  a 
%  table  of  pixel  values  through  a  pinhole  camera  model  and  return  the 
%  observed  azimuth  and  elevation  of  the  pixel  locations  in  the  image. 

%  You  will  need  to  import  a  table  of  pixel  locations  before  running. 

%%  Given: 

%  Read  in  an  image: 

original Image=imread( ’ C : \Images\matlab_images\baseline_centroids . png ’ ) ; 
%look  angles  of  Ox  Vir  converted  from  degrees  to  radians 
boresight=[214. 1906*pi/188  38 . 241944*pi/188] ; 

%  other  system  properties 
F0V=[  1.482*pi/188  1 . 857*pi/188] ; 

%%  Now  we  need  to  process  the  image  and  extract  som  values: 

%  we  require  the  size  of  the  image  in  pixels: 

%  V  H  - 

[M,  N,  d]=size(originalImage) ; 

%  Find  the  defined  variables  for  the  function  calls  below 
%Determine  the  center  points: 

0_x=(M+l)/2; 

0_y=(N+l)/2; 

%  find  virtual  focal  length  f’ 
f_prime=N/(2*atan(F0V(l , l)/2)) ; 

%%  next  create  the  "XYZ"  reference  frame  position  vector: 
d=size(streak_pixels) ; 
for  k=l :d(l , 1) 

XYZ(:,k)  =  ph_XYZ(  0_x ,0_y , streak_pixels(k , 1) , streak_pixels(k, 2) , . . . 
f_prime  ) ; 

end 

%%  Then  transform  the  XYZ  coordinate  system  to  the  SEZ  coordinate  system 
%streak  1 

SEZ  =  ph_XYZ_SEZ(  boresight(l , 1) ,pi/2-boresight(l , 2) ,XYZ  ); 

%%  Convert  the  above  LOS  vector  into  azimuth  and  elevation  spherical 
%  coordinates, 
for  k=l :d(l , 1) 

[  AZ(k, :) ,EL(k, :)]  =  az_el_third(  SEZ(:,k)  ); 

end 

%%  Detect  the  translation  present  in  the  new  image: 
translation(l , l)  =  (sum(true_location( : , 1)-AZ( : , l)))/d(l ,  1) ; 
translationCl , 2)=(sum(true_location( : ,2)-EL(: , l))3/d(l, 1) ; 

%  Translate  the  position  vectors: 
for  k=l :d(l , 1) 

[  AZEL(k,:)]  =  [AZCk, :) ,EL(k, :)]+translation; 

end 


Line  nine  is  used  to  read  in  the  given  image.  The  given  values  are  input  into  lines  11- 
13.  Line  1 1  provides  the  bore  sight  azimuth  and  elevation,  and  line  13  provides  the  FOV  of 
the  system. 
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Lines  16-24  are  used  to  process  the  given  image  to  get  required  information.  Line  18 
determines  the  number  of  pixels  in  the  x  and  y  direction.  To  continue  MATLAB  needs  to 
take  the  given  information  and  determine  the  center  point  of  the  image  {Ox,  Oy)  and  the 
virtual  focal  length  (/').  Lines  21-22  are  used  implement  Equation  2.10  to  provide  the 
center  point  of  the  image.  Line  24  implements  Equation  2. 1 1  to  determine  the  virtual  focal 
length. 

Eines  26-41  are  used  to  call  several  functions  discussed  in  Section  B.2.2.  These 
functions  loop  through  each  provided  pixel  location  and  transform  them  into  azimuth  and 
elevation  angles  as  discussed  in  Section  B.L3. 

The  section  beginning  on  line  26  uses  a  for  loop  to  create  a  position  vector  for  each 
pixel  provided  in  the  Pinhole  camera  model  reference  frame.  Eines  28-31  loop  through  the 
function  “ph_XYZ”  ,  provided  in  Section  B.2.2. 2,  to  define  the  position  vector. 

The  next  section,  beginning  on  line  33,  performs  the  transformation  to  the  SEZ 
coordinate  system  by  calling  the  function  “ph_XYZ_SEZ”  on  Eine  35.  The  function 
“phJXYZ_SEZ”  is  provided  in  Section  B. 2.2.3. 

The  final  function  call,  beginning  on  line  37,  performs  the  transformation  to  spherical 
coordinates.  Eines  39-41  loop  through  the  function  “az.eLthird”  ,  provided  in  Section 
B.2.3,  to  perform  the  transformation  for  each  position  vector  provided. 

Although  the  rotation  is  completed  by  the  “ph_XYZ_SEZ”  on  line  40,  see  Section 
B.2.5,  and  the  position  vector  translated  into  spherical  coordinates  by  line  40,  see  Section 
B.2.3,  the  position  is  still  located  at  the  origin  of  the  pinhole  camera  model.  To  translate  to 
the  correct  real  world  position  the  position  vector  requires  the  translation  given  by  Equation 
B.8.  Eines  44  and  45  determine  the  required  translation  by  averaging  the  residuals  for  both 
azimuth  and  elevation  of  the  known  stars  in  the  image.  Eines  47-49  loop  through  the 
determined  translation.  Each  angle  is  translated  as  defined  on  line  48  to  shift  to  the  correct 
position. 
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B.2.1.3  Intrinsic  Parameter  Calibration  MATLAB  Model 


The  intrinsic  parameter  calibration  model,  provided  in  Section  2.4.3,  is  implemented 
in  MATLAB  in  a  fashion  similar  to  the  pinhole  camera  model.  This  is  expected  as  they  are 
both  perspective  projection  models  [56].  The  full  script  is  available  below  in  Listing  B.3. 
Listing  B.3:  Camera  calibrated  with  intrinsic  parameter  Model  Source  Code 


%  This  script  performs  several  calculations  and  function  calls  to  run  a 
%  table  of  pixel  values  through  a  camera  with  known  intrinsic  parameters 
%  model  and  return  the  observed  azimuth  and  elevation  of  the  pixel 
%  locations  in  the  image. 

%  You  will  need  to  import  a  table  of  pixel  locations  before  running. 

%%  Given: 

%  Read  in  an  image: 

originalImage=imread( ’ C : \Images\matlab_images\baseline_centroids . png ’ ) ; 
%look  angles  of  Ox  Vir  converted  from  degrees  to  radians 
boresight=[214. 19®6*pi/180  38 . 241944*pi/180] ; 

%  other  system  properties 

F0V=[  1.402*pi/180  1 . ®57*pi/18®] ;  %  Field  of  view 
%Define  the  scaling  factors: 
s_x=l . 13 : 
s_y=®.87; 

%%  Now  we  need  to  process  the  image  and  extract  some  values: 

%  we  require  the  size  of  the  image  in  pixels: 

%  V  H  - 

[M,  N,  d]=size(originalImage) ; 

%  Find  the  defined  variables  for  the  function  calls  below 
%Determine  the  center  points: 

0_x=(M+l)/2; 

0_y=(N+l)/2; 

%  find  virtual  focal  length  f’ 
f_prime(l , l)=N/(2*atan(F0V(l , l)/2)) ; 
f_prime(l , 2)=M/(2*atanCFOV(l , 2)/2)) ; 

%%  next  create  the  "XYZ"  reference  frame  position  vector: 
d=size(streak_pixels) ; 
for  k=l :dCl , 1) 

XYZ(:,k)  =  simple_IPM(streak_pixels(k, 1) ,  streak_pixels(k, 2) ,  0_x,... 
0_y ,  s_x ,  s_y ,  f_prime  ) ; 

end 

%%  Then  transform  the  XYZ  coordinate  system  to  the  SEZ  coordinate  system 
%streak  1 

SEZ  =  XYZ_SEZ(  boresightCl, 1) ,pi/2-boresight(l,2) ,XYZ  ); 

%%  Convert  the  above  LOS  vector  into  azimuth  and  elevation  spherical 
%  coordinates, 
for  k=l :d(l , 1) 

[  AZ(k, :) ,EL(k, :)]  =  az_el_third(  SEZ(:,k)  ); 

end 

%%  Detect  the  translation  present  in  the  new  image: 
translation(l , l)  =  (sum(true_location( : , 1)-AZ(: , l)))/d(l,  1) ; 
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50  translation(l , 2)  =  (suin(true_location( :  ,2)-EL(:  ,  l)))/d(l,  1)  ; 

51 

52  %  Translate  the  position  vectors: 

53  for  k=l  :d(l ,  1) 

54  [  AZEL(k,:)]  =  [AZCk, :) ,EL(k, :)]+translation; 

55  end 


Line  10  is  used  to  read  in  the  image  that  we  are  given.  The  given  information  is 
provided  on  lines  12-17.  Line  12  provides  the  bore-sight  azimuth  and  elevation  for  the 
eenter  of  the  image.  Line  14  provides  the  FOV  covered  by  the  image.  Finally,  the  scale 
factor  is  defined  on  lines  16-17. 

Lines  20-29  are  used  to  process  the  given  image  to  get  required  information.  Line  22 
determines  the  number  of  pixels  in  the  x  and  y  direction.  To  continue  MATLAB  needs  to 
take  the  given  information  and  determine  the  center  point  of  the  image  (0^,  Oy)  and  the 
virtual  focal  length  (/')  in  both  the  vertical  and  horizontal  direction.  Lines  25-26  are  used 
implement  Equation  2.10  to  provide  the  center  point  of  the  image.  Lines  28-29  implements 
Equation  2.1 1  to  determine  the  virtual  focal  length. 

Eines  31-46  are  used  to  call  several  functions  discussed  in  Section  B.2.2.  These 
functions  loop  through  each  provided  pixel  location  and  transform  them  into  azimuth  and 
elevation  angles  as  discussed. 

The  section  beginning  on  line  3 1  uses  a  for  loop  to  create  a  position  vector  for  each 
pixel  provided  in  the  camera  model  reference  frame.  Eines  33-36  loop  through  the  function 
“simple  JPM”,  provided  in  Section  B.2.4,  to  define  the  position  vector. 

The  next  section,  beginning  on  line  38,  performs  the  transformation  to  the  SEZ 
coordinate  system  by  calling  the  function  “XYZ_SEZ”  on  Eine  40.  The  function 
“XYZ_SEZ”  is  provided  in  Section  B.2.3. 

The  next  section,  beginning  on  line  42,  performs  the  transformation  to  spherical 
coordinates.  Eines  44-46  loops  through  the  function  “az  el  third”  ,  provided  in  Section 
B.2.3,  to  perform  the  transformation. 
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Although  the  rotation  is  complete  and  the  position  vector  translated  into  spherical 
coordinates,  the  position  is  still  located  at  the  origin  of  the  calibrated  camera  model  system. 
To  translate  to  the  correct  real  world  position  the  position  vector  requires  the  translation 
given  by  Equation  B.8.  Lines  49  and  50  determine  the  required  translation  by  averaging 
the  residuals  for  both  azimuth  and  elevation  of  the  known  stars  in  the  image.  Lines  53-55 
loop  through  the  determined  translation.  Each  angle  is  translated  as  defined  on  line  54  to 
shift  to  the  correct  position  using  Equation  B.8. 

B.2.2  White  Box  Verification  of  Functions 

The  scripts  given  above  perform  several  function  calls  to  solve  parts  of  the  problem. 
Those  function  calls  are  discussed  below.  To  verify  each  function  call  it  is  compared  to  the 
theory  behind  its  creation  covered  in  Sections  2.4.1,  2.4.2,  and  2.4.3. 

B.2.2. 1  Performing  Afifine  Transform  Astrometry 

Listing  B.l  provides  all  of  the  required  information  to  perform  an  affine  transform 
on  an  image.  Line  20  of  Listing  B.l  passes  this  information  to  Listing  B.4.  Listing  B.4 
implements  the  process  provided  in  Section  2.4.1. 

Listing  B.4:  Performing  an  affine  transform 


function  [  alpha_dec  ]  =  a£fine_transform(  alpha_® , dec_8 ,  alpha,  dec,  xy,  streak  ) 
%%  Affine  Transform  Astrometry  approach  for  finding  streak  location  in  the 
%  sky 

%%  Calculate  sky  coordinates  for  reference  stars 
d=size(xy) ; 
for  k=l :d(l , 1) 

X(k, :)=Ccos(dec(k, l))*sin(alpha(k, l)-alpha_8)) . . . 

/(cos(dec_8)*cos(dec(k, l))*cos(alpha(k, l)-alpha_8) . . . 
+sin(dec_8)*sinCdec(k, 1))) ; 

Y(k, :)=C-l)*C(sin(dec_8)*cos(dec(k, l))*cos(alpha(k, l)-alpha_0) . . . 
-cos(dec_8)*sin(dec(k, 1))) . . . 

/(cos(dec_8)*cos(dec(k, l))*cos(alpha(k, l)-alpha_8) . . . 
+sin(dec_8)*sinCdec(k, 1)))) ; 

end 

%%  Solve  for  the  plate  constants 
A=zeros(d(l , 1) , 3) ; 
for  k=l :d(l , 1) 

A(k, l)=xyCk, 1) : 

A(k,2)=xy(k,2) ; 

A(k,3)=l; 

end 
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%  in  the  X  direction 
plate_x=rref ([A  X]); 

%  in  the  Y  direction 
plate_y=rref ([A  Y]); 

%%  calculate  sky  coordinates  of  each  streak 
d=size (streak) ; 

X_streak=zeros(d(l , 1) , 1) ; 

Y_streak=zeros(d(l , 1) , 1) ; 
for  k=l :dCl , 1) 

X_streak(k, :)=plate_x(l,4)*streak(k, l)+plate_x(2,4)*streak(k,2)+. . . 
plate_x(3,4) ; 

Y_streak(k, :)=plate_y(l,4)*streak(k, l)+plate_y(2,4)*streak(k,2)+. . . 
plate_y(3,4) ; 

end 

%Find  declination: 

dec_streak=asind((sin(dec_8)+Y_streak*cos(dec_8)) . . . 

./sqrt(l+X_streak.  ''2+Y_streak.  “2)) ; 

%Find  RA: 

alpha_streak=alpha_8*188/pi+atand(X_streak. . . 

. / (cos (dec_8) - Y_str eak*  sin (dec_8) ) ) ; 

alpha_dec( : , l)=alpha_streak( : , 1) ; 
alpha_dec( : , 2)=dec_streak( : , 1) ; 

end 


The  function  “affine Transform”  takes  the  inputs  inputs  alphaJ),  decS),  alpha,  dec, 
xy  and  streak  and  returns  the  look  angle  for  each  pixel  present  in  the  vector  “streak”. 
Lines  5-15  loop  through  Equations  2.2  and  2.3  to  determine  the  sky  coordinates  for  each 
of  the  reference  stars.  Lines  17-27  are  used  to  determine  the  plate  constants  by  setting  up  a 
matrix  and  finding  its  reduced  row  echelon  form,  to  solve  Equation  2.6,  providing  the  plate 
constants.  Once  the  plate  constants  are  known  Lines  33-38  loop  through  Equations  2.4  and 
2.5,  to  solve  for  the  sky  coordinates  of  each  pixel  in  the  matrix  “streak”.  Einally  Lines 
40-48  calculate  the  correct  look  angles  for  each  pixel.  Specifically  line  41  solve  Equation 
2.7  and  line  44  solves  Equation  2.8.  Lines  47  and  48  simply  place  the  returned  values  into  a 
single  matrix.  Einally,  Line  50  terminates  the  function  call  and  returns  the  matrix  of  desired 
look  angles. 

B.2.2.2  Finding  the  the  Pinhole  Camera  Model  Position  Vector 
Lines  27-31  of  Listing  B.2  loop  each  table  entry  through  the  function:  “simpleJPM” 
provided  in  Listing  B.5.  Listing  B.5  implements  Equation  2.12,  provided  in  Section  2.4.2. 
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Listing  B.5:  Determining  the  pinhole  camera  model  vector 


function  [  XYZ  ]  =  ph_XYZ(  0_x,0_y,y_prime,x_prime,f_prime  ) 

%Takes  provided  inputs  and  returns  the  pinhole  normalized  position  vector 
%  This  function  takes  inputs  x’ ,  y’ ,  0_x,  0_y,  and  f  and 
%  returns  a  pinhole  model  normalized  position  vector  using 
%  x’ ,  y’  are  the  pixel  locations  of  the  object  we  desire  to  locate  in 
%  real  space 

%  0_x,  Y_x  are  the  center  pixel  coordinates  of  the  image. 

%  f’  is  the  focal  length  of  the  camera  system 
X=0_x-x_prime ; 

Y=y_prime-0_y ; 

Z=f_prime ; 

XYZ=[X;  Y;  Z] ; 
end 


The  function  “phJXYZ”  takes  the  inputs  inputs  x' ,  y'.  Ox,  Oy,  and  f  and  returns  the 
pinhole  normalized  position  vector.  Lines  9-11  solve  for  the  first,  second,  and  third  entry 
of  Equation  2.12  respectively.  Line  13  terminates  the  function  call  and  returns  the  \XYZY 
vector. 

B.2.2.3  Rotating  the  Pinhole  Camera  Model  Vector 
The  section  beginning  on  line  33  of  Listing  B.2,  implements  the  rotation  discussed  in 
Section  B.l.l.  Line  35  calls  the  function  “phJXYZ  SEZ”,  provided  in  Eisting  B.6  to  rotate 
the  position  vectors  to  align  it  with  the  SEZ  vector. 

Eisting  B.6:  Rotating  the  camera  model  vector  to  SEZ 


function  [  SEZ  ]  =  ph_XYZ_SEZC  theta, phi , XYZ  ) 

%Rotates  a  normalized  position  vector  to  align  it  with  SEZ 
%  This  function  takes  inputs  phi,  theta,  and  the  vector  XYZ  and  returns 

%  the  SEZ  vector  after  performing  a  2-3  rotation  on  it. 

%  phi  is  the  angle  of  the  first  rotation  in  radians 
%  theta  is  the  angle  of  the  second  rotation  in  radians 
R0T2  =  [cos(phi)  0  -sin(phi) ; . . . 

010;  sin(phi)  0  cos(phi)]: 

ROTS  =  [cos (theta)  sinCtheta)  0;... 

-sin(theta)  cos(theta)  0;  001]; 

SEZ  =  ROTS*  R0T2*  XYZ; 
end 


The  function  “ph_XYZ_SEZ”  takes  inputs  theta,  phi,  and  the  vector  \XYZY  and 
returns  the  SEZ  vector  after  performing  a  2-3  rotation  on  it.  Eines  7-8  defines  the  rotation 
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given  by  Equation  B.l,  and  lines  9-10  define  the  rotation  given  by  Equation  B.2.  Einally 
line  11  performs  the  2-3  rotation,  given  by  Equation  B.3,  on  the  provided  veetor  and  line 
12  terminates  the  funetion  eall  and  returns  the  [S EZ\^  veetor. 

B.2.3  Transforming  into  Spherical  Coordinates 

Both  Eisting  B.2  and  Eisting  B.3,  perform  the  transformation  to  spherieal  eoordinates. 
This  is  aeeomplished  through  the  funetion  “az  el  third”  to  perform  the  transformation. 

The  funetion  “az_el_third”,  provided  in  Eisting  B.7,  takes  the  provided  SEZ  veetor 
and  returns  the  spherieal  eoordinates  (J3,  el).  Eine  9  implements  Equation  B.6.  Eine  10 
implements  Equation  B.7.  Einally  line  1 1  terminates  the  funetion  eall  and  returns  the  f  and 
el  angles. 

Eisting  B.7:  Transforming  SEZ  to  look  angle 


function  [  AZ,  EL  ]  =  az_el_third(  rho  ) 

%Takes  provided  inputs  and  returns  the  spherical  coordinates  (AZ,  EL) 
%  This  function  takes  input  vector  rho  and  returns  the  azimuth  and 
%  elevation  associated  with  the  provided  vector. 

%  rho  is  a  position  vector  in  the  SEZ  coordinate  system. 

%  Az  is  the  azimuth  of  the  angle  returned  in  degrees 
%  el  is  the  elevation  of  the  angle  returned  in  degrees 

AZ=atand(rho(2 , l)/rho(l , 1))+18® ;  %(-(-))/+  is  third  quad  +18® 

EL=asind(-rhoC3 , l)/normCrho)) ; 

end 


B.2.4  Finding  the  Calibrated  Camera  Model  Position  Vector 
Eines  33-36  of  Eisting  B.3,  loop  through  a  function,  “simple  JPM”  provided  in  Eisting 
B.8.  Eisting  B.8  implements  Equation  2.14,  provided  in  Seetion  2.4.3. 

Eisting  B.8:  Determining  the  ealibrated  eamera  veetor 


function  [  XYZ  ]  =  simple_IPMC  x_prime,  y_prime,  0_x,  0_y,  s_x,  s_y,  f  ) 
%Takes  provided  inputs  and  returns  the  simplified  answer  position  vector 
%  This  function  takes  inputs  x’ ,  y’ ,  0_x,  0_y,  s_x,  s_y,  and  f  and 

%  returns  a  skew  free  normalized  position  vector  using  a  camera  with 

%  intrensic  parameters  model. 

%  x’ ,  y’  are  the  pixel  locations  of  the  object  we  desire  to  locate  in 
%  real  space 

%  0_x,  Y_x  are  the  center  pixel  coordinates  of  the  image. 

%  s_x,  s_y  are  the  scale  factors  for  the  image 
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%  f  is  the  focal  length  of  the  camera  system 


X=Cx_prime-0_x)/(s_x*f) ; 
Y=Cy_prime-0_y)/(s_y*f) ; 
Z=l; 

XYZ=[X;  Y;  Z]  ; 
end 


The  function  “simpleJPM”  takes  the  inputs  inputs  x' ,  y' ,  Ox,  Oy,  Sx,  Sy,  and  /  and 
returns  a  skew  free  normalized  position  vector  using  a  camera  with  intrinsic  parameters 
model.  Line  one  initializes  the  function.  Lines  2-10  provide  a  description  of  the  function 
and  how  it  works.  Lines  12,  13,  and  14  solve  for  the  first,  second,  and  third  entry  of 
Equation  2.14  respectively.  Line  18  terminates  the  function  call  and  returns  the  \XYZY 
vector. 

B.2.5  Rotating  the  Calibrated  Camera  Model  Vector 

Listing  B.3  calls  the  function  “XYZ_SEZ”  to  implement  the  rotation  discussed  in 
Section  B.  1.2  to  align  the  position  vector  with  the  SEZ  vector.  The  function  “XYZ_SEZ”, 
provided  in  Listing  B.9,  takes  inputs  phi,  theta,  and  the  vector  XYZ  and  returns  the  SEZ 
vector  after  performing  a  1-3  rotation  on  it.  Lines  9-10  define  the  rotation  given  by  Equation 
B.4,  and  lines  11-12  define  the  rotation  given  by  Equation  B.2.  Einally  line  14  performs 
the  1-3  rotation,  given  by  Equation  B.5,  on  the  provided  vector  and  line  16  terminates  the 
function  call  and  returns  the  \SEZY  vector. 

Listing  B.9:  Rotating  the  calibrated  camera  vector  to  SEZ 


function  [  SEZ  ]  =  XYZ_SEZ(  phi , theta, XYZ  ) 

%Rotates  a  normalized  position  vector  to  align  it  with  SEZ 
%  This  function  takes  inputs  phi,  theta,  and  the  vector  XYZ  and  returns 

%  the  SEZ  vector  after  performing  a  1-3  rotation  on  it. 

%  phi  is  the  angle  of  the  first  rotation  in  radians 
%  theta  is  the  angle  of  the  second  rotation  in  radians 

%Convert  Normalized  Coordinates  to  SEZ  coordinates: 

ROTl  =  [1  ®  Q;  0  cos(phi)  sinCphi) ; . . . 

0  -sin(phi)  cos(phi)] ; 

ROTS  =  [cos (theta)  sin(theta)  0;... 

-sin(theta)  cos(theta)  0;  001]; 
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14  SEZ  =  R0T3*  ROTl*  [XYZ]  ; 

15 

16  end 


B.2.5.1  20  Star  Verification 

The  Second  Verification  process  compares  the  expected  look  angles  of  20  stars  with 
known  look  angles.  This  was  accomplished  by  generating  a  noise  free  simulated  image 
of  a  small  FOV  of  the  night  sky.  The  image  was  processed  to  extract  the  pixel  locations 
of  the  known  stars,  and  those  pixel  locations  were  used  as  inputs  to  each  astrometry  tool. 
The  results,  provided  graphically  can  be  seen  in  Figures  B.9,  B.IO,  and  B.ll  are  directly 
compared  to  the  true  look  angle  of  each  star.  Each  model  is  discussed  in  more  detail  below. 

The  results  of  the  affine  transform  verification  are  shown  in  Figure  B.9,  again  directly 
overlying  the  true  positions  for  visual  comparison.  As  can  be  seen  in  Figure  B.9  the  overlay 
of  the  true  position  on  the  returned  values  is  almost  an  exact  match.  As  the  affine  transform 
method  attempts  to  account  for  all  distortion  this  outcome  is  as  expected.  Figure  B.9  further 
verifies  that  the  MATLAB  script  provided  in  Listing  B.l,  is  functioning  as  expected. 
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Figure  B.9:  Affine  transform  astrometry  overlying  true  positions 


The  results  of  the  pinhole  camera  verification  are  shown  in  Figure  B.IO,  again 
overlying  the  true  positions,  for  visual  comparison.  As  shown  in  Figure  B.IO,  the  pinhole 
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camera  results  are  best  near  the  eenter  of  the  image  and  deviate,  mostly  in  elevation,  as  we 
deviate  from  eenter.  This  focal  length  distortion  is  not  unexpeeted  as  the  pinhole  model  only 
aeeounts  for  the  focal  length  in  the  horizontal  direetion,  and  unlike  a  model  that  ineludes 
intrinsie  parameters,  does  not  aeeount  for  seale  in  the  image.  Seale  is  present  even  in  the 
noise  free  image  as  the  CCD  pixels  are  reetangular  and  not  perfeetly  square.  As  expeeted 
the  further  from  the  prineiple  point  we  get  the  more  the  seale  affects  the  image.  As  the 
distortions  in  the  returned  positions  ean  be  aceounted  for  in  the  model.  Figure  B.IO  further 
verifies  that  the  MATLAB  script  provided  in  Listing  B.2  is  functioning  correctly. 
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Figure  B.IO:  Pinhole  camera  model  overlying  true  positions 


Finally,  the  results  of  the  intrinsie  parameter  calibrated  eamera  model  verifieation  are 
shown  in  Figure  B.ll,  again  overlying  the  true  positions,  for  visual  eomparison.  For  the 
intrinsie  eamera  ealibrated  image  the  overlay  of  the  true  position  on  the  returned  values  is 
almost  an  exaet  mateh.  As  the  ealibrated  image  aeeounts  for  all  of  the  eamera  distortion 
you  would  expect  this  result  to  be  true  for  a  noise  free  image.  Figure  B.ll  further  verifies 
that  the  MATLAB  script  provided  in  Listing  B.3,  is  funetioning  as  expected. 
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Figure  B.ll:  Intrinsic  parameters  calibrated  model  overlying  true  positions 


B.2.6  Validation  of  the  Angle  Extraction  Tools 

To  validate  each  of  the  astrometry  method  scripted  in  Section  B.2.1  requires 
confirmation  that  when  the  required  inputs  are  given,  look  angles  for  the  desired  pixel 
locations  are  returned.  This  was  demonstrated  in  Section  B.2.5.1  as  each  method  was 
verified.  Therefore,  further  validation  of  the  models  is  not  required. 
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Appendix  C:  Orbital  Determination  Unscented  Kalman  Filter 


The  focus  of  this  dissertation  was  on  the  extraction  of  angles  data  from  images 
provided  by  external  sensors.  The  method  developed  isolated  streaks  and  extracted  angles 
data  as  described  in  Appendix  A.  The  intention  was  to  be  able  to  pass  the  angles  data  to 
an  existing  orbital  algorithm  including:  Laplace,  Gauss,  and  the  Double-R  method  using 
standard  orbital  determination  algorithms  [8]. 

When  the  determination  of  orbital  data  was  required  for  Phase  II  Gary  Goff’s  UKF 
algorithm  was  used.  This  appendix  describes  how  that  algorithm  was  implemented  for 
the  work  presented  here.  The  majority  of  this  Appendix  was  written  by  Gary  Goff  as  his 
contribution  to  the  Paper:  Updating  Position  Data  from  Full  Serendipitous  Satellite  Streaks. 


C.l  Orbital  Determination  UKF 

When  performing  an  orbit  estimation,  a  key  early  decision  is  to  determine  an  element 
set  to  represent  both  the  orbital  state  and  the  covariance.  An  orbital  solution  is  composed 
of  n  =  6  components,  where  n  is  the  size  of  the  state  vector.  A  traditional  choice  is 
Earth-centered  Cartesian  coordinates  of  position  and  velocity,  as  the  covariance  is  easy 
to  visualize.  Recent  works  by  the  team  at  Numerica,  Inc.  have  shown  in  several  papers 
that  equinoctial  elements  better  preserve  the  underlying  error  distribution  [73,  94,  95]. 
Based  on  these  works  and  initial  testing  by  the  authors,  an  equinoctial  element  UKF 
was  implemented  to  process  telescope  observations  while  estimating  the  orbital  state  and 
covariance.  The  state  is  define  as  a  set  of  equinoctial  elements  [8]. 

r  ^T 


X  = 


a  h  k  pe  qe  £ 


(C.l) 


Unfortunately,  the  elements  are  not  easily  visualized  like  classical  orbital  elements, 
but  equinoctial  elements  are  often  used  to  avoid  singularities  present  in  classical  orbital 
elements  [96] .  The  estimated  covariance  matrix  for  the  corresponding  state  is  captured  in 
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the  matrix  P.  The  UKF  routine  requires  an  initial  state  and  covariance  which  is  provided 
from  a  previous  estimate  or  determined  via  initial  orbit  determination  routines  for  newly 
discovered  RSOs. 


As  mentioned  in  the  background  section,  the  UKF  implements  the  UT  which  uses 
weighted  sampled  sigma  points  (specifically  2n+ 1)  to  estimate  the  nonlinear  transformation 
[97].  The  symmetric  sigma  point  weighting  approach  is  commonly  used  and  the  optimal 
K  value  for  Gaussian  distributions  is  defined  such  that  n  +  k  =  3  [98].  The  mean  and 
covariance  weights  of  sigma  point  j  are  defined  as: 

K 


= 

Wyn  C 


n  K 


1 


for  j  =  l,..,2n 


(C.2) 

(C.3) 


2{n  +  k) 

where  a:  =  -3  is  recommended  for  orbital  elements. 

The  UT  requires  taking  the  square  root  of  the  covariance  matrix.  For  a  symmetric, 
positive  definite  matrix  (true  for  orbital  covariance  matrices),  this  operation  is  easily 
performed  via  a  Cholesky  decomposition  [99].  Implementing  the  UT  into  the  UKF  for 
orbit  estimation  requires  two  main  steps:  propagate  and  update.  The  detailed  steps  of  the 
orbital  estimation  UKF  are  provided  in  Algorithm  1  [97]. 

The  UKF  in  Algorithm  1  requires  propagating  each  state  sigma  point  using  some 
form  of  orbital  state  propagation  (General  or  Special  Perturbation).  For  this  research. 
Microcosm’s  HPOP  propagator  is  implemented  due  to  its  speed  and  accuracy  [100].  Once 
all  sigma  points  are  propagated,  the  propagated  covariance,  P,  is  estimated  within  the 
UKF.  The  process  noise  covariance,  Q,  captures  the  distribution  of  propagation  errors. 
Due  to  relatively  short  propagation  times  and  orbit  fitting  test  results,  Q  was  set  to  zero 
as  the  propagation  step  did  not  induce  significant  errors.  The  actual  observation  of  the 
state  is  performed  at  the  telescope  when  it  reports  a  right  ascension,  a,  and  declination,  6, 


angle.  The  actual  observation  vector  is  defined  such  that  y  = 


a  6 


.  The  telescope’s 


errors  are  assumed  to  be  an  uncorrelated  Gaussian  with  a  covariance  captured  in  R  such 
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that  R  =  diag  (cTq,,  cr^),  where  cr  represents  the  standard  deviation  of  the  error.  The  bias 
error  is  discussed  separately  in  the  calibration  section  of  this  paper.  Once  all  sigma  points 
are  propagated,  each  is  transformed  into  a  predicted  observation  through  the  function  G. 
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Step  7  in  Algorithm  1,  the  propagated  sigma  point  states  are  transformed  from  equinoctial 
elements  into  an  anticipated  observation  vector  y  by  leveraging  Analytical  Graphics,  Inc.’s 
System  Tool  Kit  ®.  The  residual,  v,  is  defined  as  the  difference  between  the  observed 
and  the  calculated  observation.  Since  the  posterior  probability  distribution  function  (pdf) 
covariance  update,  P,  within  a  Kalman  filter  is  independent  of  the  residual,  it  is  important  to 
assess  the  relationship  between  the  covariance  and  residual.  Since  the  matrix  S  represents 
the  predicted  observation  covariance,  'T  defines  the  squared  Mahalanobis  distance  of  the 
observation  and  indicates  the  accuracy  of  the  orbital  fit  [16,  101,  102].  The  remaining 
equations  in  Algorithm  1  follow  a  traditional  Kalman  filter  and  result  in  an  updated  estimate 
for  the  state,  x,  and  covariance,  P. 

C.2  Starting  Covariance 

To  implement  the  UKF,  a  starting  covariance  is  needed  to  accompany  the  TLE  of 
interest.  Since  the  candidate  RSO  is  in  a  GEO  orbit  and  observed  frequently,  a  high 
accuracy  covariance  as  proposed  by  Horwood  is  used  [73].  Define  the  starting  covariance 
in  equinoctial  elements  such  that  Pq  =  A^A,  where 

A  =  diag(0.05  km,  10-^  10-^  10-^  10■^97  x  10“^  rad)  (C.4) 

The  initial  state  is  pulled  directly  from  the  published  TEE. 

C.3  Data  Sets 

Eor  the  analysis  in  Chapter  4,  the  comparison  of  orbital  data  generated  using  three 
subsets  of  the  data  was  required.  The  first  set  of  data  is  created  strictly  from  the  control 
group,  to  provide  a  baseline  analysis.  The  second  data  set  generates  an  orbit  by  excluding 
data  from  the  night  with  the  simulated  unbounded  streak.  The  second  data  set  provide 
estimated  orbital  positions  assuming  that  the  captured  unbounded  streak  was  thrown  out, 
and  the  orbit  was  determined  using  only  the  two  nights  of  data.  The  final  data  set 
incorporates  the  unbounded  streak  data.  These  data  sets  allow  for  analysis  of  the  usefulness 
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of  the  position  data  in  the  unbounded  streak.  The  impact  of  using  the  unbounded  streak  vs. 
ignoring  the  streak  and  using  fewer  observations  is  discussed  below. 
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Appendix  D:  Useful  Tools 


D.l  Software 

D.Ll  MATLAB 

The  chief  advantage  of  MATLAB  is  its  ability  to  execute  high-level  matrix  algebra 
and  built  in  data  processing  routines  [103].  This  benefit  makes  MATLAB  a  very  capable 
image  processing  tool.  It  is  capable  of  performing  complex  matrix  calculations  as  well  as 
output  plots,  images,  and  tables  of  data  [103]. 

The  built  in  image  processing  tool  suite  is  capable  of  analyzing  fits  files  and  extracting 
information  such  as  light  intensity,  pixel  location,  Gaussian  spread  functions,  etc.  Currently 
this  technique  is  the  most  promising  image  processing  tool  used  for  this  research  and  will 
continue  to  be  used  in  the  absence  of  better  alternatives. 

Vallado  also  translated  his  astrodynamics  algorithms  infrom  fortran  to  MATLAB 
[8].  As  this  research  hopes  to  demonstrate  the  applicability  of  processing  external  data 
for  astrodynamics  application  Vallado’s  work  may  prove  a  suitable  tool  for  the  orbital 
determination  and  updates. 

D.1.2  Smithsonian  Astrophysical  Observatorylmage  ds9 

Smithsonian  Astrophysical  Observatorylmage  ds9  (ds9)  is  a  powerful,  open  source 
tool  available  for  Windows,  Linux,  and  Macs  developed  by  astronomers  at  Harvard.  ds9 
allows  scientists  to  make  quantitative  and  valid  measurements  from  fits  files.  It  can  handle 
every  type  of  fits  files  but  prefers  images;  if  used  for  spectra  most  astronomers  will  use 
other  tools  in  conjunction.  ds9  is  also  widely  used  by  astronomers  worldwide  to  manipulate 
images  [22]. 

Image  World  Coordinate  System  Setting  Utility  (imwcs)  is  a  program  written  by 
Smithsonian  Astrophysical  Observatory  (SAO)  to  work  within  ds9  to  automatically  overlay 
a  grid  of  a  and  6  onto  a  fits  image.  This  technique  is  accomplished  by  embedding  the  WCS 
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from  a  star  catalog  into  the  fits  file  header.  After  the  WCS  has  been  written  to  the  file 
any  fits  processing  tool  can  be  used  to  process  the  image  and  extract  the  a  and  S  for  the 
satellite  streak.  Unfortunately  this  program  only  runs  in  an  Image  Reduction  and  Analysis 
Facility  (IRAF)  environment  limiting  the  types  of  computers  that  can  use  it. 

D.I.3  Teletrak 

AFIT’s  teletrak  system  was  originally  conceived  to  determine  the  orbit  of  LEO 
satellites  using  commercial  telescopes.  It  contains  a  MATLAB  Graphical  User  Interface 
(GUI)  interface.  Once  setup  the  GUI  displays  a  comprehensive  star  map  for  the  location 
of  the  telescope  used  along  with  brightness  predictions  and  satellite  targeting  information 
[11]. 

This  targeting  information  can  be  used  to  predict  generated  satellite  streaks  within  the 
telescope’s  FOV.  The  only  potential  drawback  is  that  some  modification  would  be  required 
to  get  the  GUI  to  behave  as  any  telescope  not  already  included  in  the  teletrak  network  [11]. 

D.1.4  Systems  Tool  Kit 

Analytical  Graphics  Inc.  (AGI)’s  STK  suite  is  designed  to  simulate  space  systems. 
Ground  sites  can  be  created,  including  sensor  FOV,  and  contact  reports  generated  for  all 
satellites  with  public  TLE’s,  reporting  azimuth  and  elevation  for  all  times  the  vehicle  has 
line  of  sight.  Combined  with  the  ability  to  calculate  sun  angles,  a  script  can  be  written  to 
effectively  simulate  satellite  streaks  across  a  telescope  [104]. 

D.1.5  Stellarium  and  KStars 

Stellarium  and  KStars  are  software  suites  designed  to  allow  a  user  to  use  a  home 
Personal  Computer  (PC)  as  a  virtual  planetarium.  They  calculate  the  position  of 
astronomical  objects  and  events  and  simulate  the  sky  from  a  user  defined  location.  They 
can  serve  as  both  an  educational  tool  as  well  as  observational  aids  for  astronomers.  As  an 
observational  aid  either  KStars  or  Stellarium  can  be  setup  to  imitate  a  telescope  system  as 
long  as  the  user  provides:  resolution  x  and  y  pixel  count,  the  chip  width  and  height  in  mm. 
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and  the  native  focal  length.  Along  with  the  value  of  the  pixel  width  and  height  in  microns, 
the  multiplier  of  any  lens  used,  and  the  diameter  of  the  telescope  used  [86,  105]. 

Some  of  the  advantages  of  Stellarium  include  automated  access  to  outside  sources, 
such  as  celestrack,  to  upload  the  publicly  available  TLE.  Using  this  TLE  it  can  determine 
if  a  satellite  is  within  the  EOV  of  the  telescope,  the  sun  angle,  and  its  visual  magnitude  to 
determine,  and  automatically  report,  if  the  RSO  is  visible  to  the  observer  [86]. 

Stellarium  also  contains  an  Application  Programming  Interface  (API)  scripting 
interface  allowing  it  to  be  plugged  into  a  larger  modeling  tool,  such  as  Phoenix  Integration’s 
Model  Center,  to  simulate  a  single  node  in  a  larger  sensor  network.  This  scripting  language 
can  also  be  used  with  the  above  capabilities  to  model  what  a  satellite  streak  on  a  CCD  by 
scripting  the  sensor  characteristics  and  exposure  time  to  report  any  visible  satellites  a  and 
5  in  that  time  frame  [86]. 

Conversely  KStars  allows  for  the  creation  of  artificial  satellites,  similar  to  Stellarium, 
using  TEE  data  but  it  is  unclear  from  the  documentation  if  it  includes  a  sun  angle  calculation 
to  determine  visibility  through  the  telescope  [105]. 

D.2  Hardware 

D.2.1  Meade  Instruments 

Several  of  the  telescopes  used  at  AEIT  are  manufactured  by  Meade  Instrument 
Corporation,  including  both  a  10  inch  and  16  inch  telescope,  currently  mounted  on  Meade 
Autostar  II  mounts.  These  mounts  allow  an  operator  to  align  the  optical  system  using  one 
or  two  stars  contained  within  the  internal  catalog,  point  the  telescope  at  a  desired  point  in 
the  sky,  and  track  an  object  [106].  Past  work  at  AEIT  has  created  programs  to  adjust  the 
tracking  rate  of  the  scope  so  that  it  will  not  only  track  at  the  rate  of  the  celestial  sphere, 
useful  for  astrometric  observations,  but  can  also  track  a  EEO  satellite  as  it  passes  overhead 
[11]. 
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